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ESSENTIAL  ELEMENTS  OF  COMPUTATIONAL  ALGORITHMS  FOR  AERODYNAMIC 

ANALYSIS  AND  DESIGN* 

ANTONY  JAMESON^ 


Abstract.  This  paper  traces  the  development  of  computational  fluid  dynamics  as  a  tool  for  aircraft 
design.  It  addresses  the  requirements  for  effective  industrial  use,  and  trade-offs  between  modeling  accuracy 
and  computational  costs.  Essential  elements  of  algorithm  design  are  discussed  in  detail,  together  with  a 
unified  approach  to  the  design  of  shock  capturing  schemes.  Finally,  the  paper  discusses  the  use  of  techniques 
drawn  from  control  theory  to  determine  optimal  aerodynamic  shapes.  In  the  future  multidisciplinary  analysis 
and  optimization  should  be  combined  to  provide  an  integrated  design  environment. 
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Subject  classification.  Fluid  Mechanics 

1.  Introduction.  Computational  methods  first  began  to  have  a  significant  impact  on  aerodynamics 
analysis  and  design  in  the  period  of  1965-75.  This  decade  saw  the  introduction  of  panel  methods  which 
could  solve  the  linear  flow  models  for  arbitrarily  complex  geometry  in  both  subsonic  and  supersonic  flow 
[68, 164,  200].  It  also  saw  the  appearance  of  the  first  satisfactory  methods  for  treating  the  nonlinear  equations 
of  transonic  flow  [135,  134,  75,  76,  51,  63],  and  the  development  of  the  hodograph  method  for  the  design  of 
shock  free  supercritical  airfoils  [19]. 

Computational  Fluid  Dynamics  (CFD)  has  now  matured  to  the  point  at  which  it  is  widely  accepted  as  a 
key  tool  for  aerodynamic  design.  Algorithms  have  been  the  subject  of  intensive  development  for  the  past  two 
decades.  The  principles  imder lying  the  design  and  implementation  of  robust  schemes  which  can  acciurately 
resolve  shock  waves  and  contact  discontinuities  in  compressible  flows  are  now  quite  well  established.  It  is 
also  quite  well  understood  how  to  design  high  order  schemes  for  viscous  flow,  including  compact  schemes 
and  spectral  methods.  Adaptive  refinement  of  the  mesh  interval  (h)  and  the  order  of  approximations  (p) 
has  been  successfully  exploited  both  separately  and  in  combination  in  the  h-p  method  [138].  A  continuing 
obstacle  to  the  treatment  of  configurations  with  complex  geometry  has  been  the  problem  of  mesh  generation. 
Several  general  techniques  have  been  developed,  including  algebraic  transformations  and  methods  based  on 
the  solution  of  elHptic  and  hyperbohc  equations.  In  the  last  few  years  methods  using  unstructured  meshes 
have  also  begun  to  gain  more  general  acceptance.  The  Dassault-INRIA  group  led  the  way  in  developing  a 
finite  element  method  for  transonic  potential  flow.  They  obtained  a  solution  for  a  complete  Falcon  50  as 
early  as  1982  [29].  Euler  methods  for  unstructured  meshes  have  been  the  subject  of  intensive  development 
by  several  groups  since  1985  [118,  90,  89,  183,  18],  and  Navier-Stokes  methods  on  unstructured  meshes  have 
also  been  demonstrated  [126,  127,  17]. 

Despite  the  advances  that  have  been  made,  CFD  is  still  not  being  exploited  as  effectively  as  one  would 
hke  in  the  design  process.  This  is  partly  due  to  the  long  set-up  and  high  costs,  both  human  and  computational 
of  complex  flow  simulations.  The  essential  requirements  for  industrial  use  are: 

1.  assured  accuracy 
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2.  acceptable  computational  and  human  costs 

3.  fast  turn  around. 

Improvements  are  still  needed  in  all  three  areas.  In  particular,  the  fidelity  of  modelling  of  high  Reynolds 
number  viscous  flows  continues  to  be  limited  by  computational  costs.  Consequently  accurate  and  cost 
effective  simulation  of  viscous  flow  at  Reynolds  numbers  associated  with  full  scale  flight,  such  as  the  prediction 
of  high  lift  devices,  remains  a  challenge.  Several  routes  are  available  toward  the  reduction  of  computational 
costs,  including  the  reduction  of  mesh  requirements  by  the  use  of  higher  order  schemes,  improved  convergence 
to  a  steady  state  by  sophisticated  acceleration  methods,  fast  inversion  methods  for  implicit  schemes,  and  the 
exploitation  of  massively  parallel  computers. 

Another  factor  limiting  the  effective  use  of  CFD  is  the  lack  of  good  interfaces  to  computer  aided  design 
(CAD)  systems.  The  geometry  models  provided  by  existing  CAD  systems  often  fail  to  meet  the  requirements 
of  continuity  and  smoothness  needed  for  flow  simulation,  with  the  consequence  that  they  must  be  modified 
before  they  can  be  used  to  provide  the  inpxit  for  mesh  generation.  This  bottleneck,  which  impedes  the 
automation  of  the  mesh  generation  process,  needs  to  be  ehminated,  and  the  CFD  software  should  be  fully 
integrated  in  a  numerical  design  environment.  In  addition  to  more  accurate  and  cost-effective  flow  prediction 
methods,  better  optimizations  methods  are  also  needed,  so  that  not  only  can  designs  be  rapidly  evaluated, 
but  directions  of  improvement  can  be  identified.  Possession  of  techniques  which  result  in  a  faster  design 
cycle  gives  a  crucial  advantage  in  a  competitive  environment. 

A  critical  issue,  examined  in  the  next  section,  is  the  choice  of  mathematical  models.  What  level  of 
complexity  is  needed  to  provide  sufficient  accuracy  for  aerodjmamic  design,  and  what  is  the  impact  on  cost 
and  turn-around  time?  Section  3  addresses  the  design  of  numerical  algorithms  for  flow  simulation.  Section 
4  presents  the  results  of  some  numerical  calculations  which  require  moderate  computer  resources  and  could 
be  completed  with  the  fast  turn-around  required  by  industrial  users.  Section  5  discusses  automatic  design 
procedmes  which  can  be  used  to  produce  optimum  aerod3mamic  designs.  Finally,  Section  6  offers  an  outlook 
for  the  future. 

2.  The  Complexity  of  Fluid  Flow  and  Mathematical  Modeling. 

2.1.  The  Hierarchy  of  Mathematical  Models.  Many  critical  phenomena  of  fluid  flow,  such  as  shock 
waves  and  turbulence,  are  essentially  non-linear.  They  also  exhibit  extreme  disparities  of  scales.  While  the 
actual  thickness  of  a  shock  wave  is  of  the  order  of  a  mean  free  path  of  the  gas  particles,  on  a  macroscopic 
scale  its  thickness  is  essentially  zero.  In  turbulent  flow  energy  is  transferred  from  large  scale  motions  to 
progressively  smaller  eddies  until  the  scale  becomes  so  small  that  the  motion  is  dissipated  by  viscosity.  The 
ratio  of  the  length  scale  of  the  global  flow  to  that  of  the  smallest  persisting  eddies  is  of  the  order  Re  ^ ,  where 
Re  is  the  Reynolds  number,  typically  in  the  range  of  30  million  for  an  aircraft.  In  order  to  resolve  such  scales 
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in  all  three  space  directions  a  computational  grid  with  the  order  of  Re^  cells  would  be  required.  This  is 
beyond  the  range  of  any  current  or  foreseeable  computer.  Consequently  mathematical  models  with  varying 
degrees  of  simplification  have  to  be  introduced  in  order  to  make  computational  simulation  of  flow  feasible 
and  produce  viable  and  cost-effective  methods. 

Figure  1  (supplied  by  Pradeep  Raj)  indicates  a  hierarchy  of  models  at  different  levels  of  simplification 
which  have  proved  useful  in  practice.  Efficient  flight  is  generally  achieved  by  the  use  of  smooth  and  stream¬ 
lined  shapes  which  avoid  flow  separation  and  minimize  viscous  effects,  with  the  consequence  that  useful 
predictions  can  be  made  using  inviscid  models.  Inviscid  calculations  with  boundary  layer  corrections  can 
provide  quite  accurate  predictions  of  lift  and  drag  when  the  flow  remains  attached,  but  iteration  between  the 
inviscid  outer  solution  and  the  inner  boundary  layer  solution  becomes  increasingly  difficult  with  the  onset 
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of  separation.  Procedures  for  solving  the  full  viscous  equations  are  likely  to  be  needed  for  the  simulation  of 
arbitrary  complex  separated  flows,  which  may  occur  at  high  angles  of  attack  or  with  bluff  bodies.  In  order 
to  treat  flows  at  high  Reynolds  numbers,  one  is  generally  forced  to  estimate  turbulent  effects  by  Re3niolds 
averaging  of  the  fluctuating  components.  This  requires  the  introduction  of  a  turbulence  model.  As  the 
available  computing  power  increases  one  may  also  aspire  to  large  eddy  simulation  (LES)  in  which  the  larger 
scale  eddies  are  directly  calculated,  while  the  influence  of  turbulence  at  scales  smaller  than  the  mesh  interval 
is  represented  by  a  subgrid  scale  model. 


Fig.  1.  Hierarchy  of  Fluid  Flow  Models 


2.2.  Computational  Costs.  Computational  costs  vary  drastically  with  the  choice  of  mathematical 
model.  Panel  methods  can  be  effectively  used  to  solve  the  linear  potential  flow  equation  with  higher-end 
personal  computers  (with  an  Intel  Pentium  microprocessor,  for  example).  Studies  of  the  dependency  of  the 
result  on  mesh  refinement,  performed  by  this  author  and  others,  have  demonstrated  that  inviscid  transonic 
potential  flow  or  Euler  solutions  for  an  airfoil  can  be  accurately  calculated  on  a  mesh  with  160  cells  around 
the  section,  and  32  cells  normal  to  the  section.  Using  multigrid  techniques  10  to  25  cycles  are  enough  to 
obtain  a  converged  result.  Consequently  airfoil  calculations  can  be  performed  in  seconds  on  a  Cray  YMP,  and 
can  also  be  performed  on  Pentium-class  personal  computers.  Correspondingly  accurate  three-dimensional 
inviscid  calculations  can  be  performed  for  a  wing  on  a  mesh,  say  with  192x32x48  =  294, 912  cells,  in  about 
5  minutes  on  a  single  processor  Cray  YMP,  or  less  than  a  minute  with  eight  processors,  or  in  less  than  an 
hour  on  a  workstation  such  as  a  Sihcon  Graphics  Indigo  2. 

Viscous  simulations  at  high  Reynolds  numbers  require  vastly  greater  resources.  Careful  two-dimensional 
studies  of  mesh  requirements  have  been  carried  out  at  Princeton  by  Martinelli  [122] .  He  found  that  on  the 
order  of  32  mesh  intervals  were  needed  to  resolve  a  turbulent  boundary  layer,  in  addition  to  32  intervals  be¬ 
tween  the  boundary  layer  and  the  far  field,  leading  to  a  total  of  64  intervals.  In  order  to  prevent  degradations 
in  accuracy  and  convergence  due  to  excessively  large  aspect  ratios  (in  excess  of  1,000)  in  the  surface  mesh 
cells,  the  chordwise  resolution  must  also  be  increased  to  512  intervals.  Reasonably  accurate  solutions  can  be 
obtained  in  a  512x64  mesh  in  100  multigrid  cycles.  TVanslated  to  three  dimensions,  this  would  imply  the 
need  for  meshes  with  5-10  million  cells  (for  example,  512x64x256  =  8,388,608  cells  as  shown  in  Figure  2). 
When  simulations  are  performed  on  less  fine  meshes  with,  say,  500,000  to  1  million  cells,  it  is  very  hard  to 
avoid  mesh  dependency  in  the  solutions  as  well  as  sensitivity  to  the  turbulence  model. 

A  typical  algorithm  requires  of  the  order  of  5,000  floating  point  operations  per  mesh  point  in  one 
multigrid  iteration.  With  10  million  mesh  points,  the  operation  count  is  of  the  order  of  0.5x10^^  per  cycle. 
Given  a  computer  capable  of  sustaining  10^^  operations  per  second  (100  gigaflops),  200  cycles  could  then  be 
performed  in  100  seconds.  Simulations  of  unsteady  viscous  flows  (flutter,  buffet)  would  be  likely  to  require 
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Fig.  2.  Mesh  Requirements  for  a  Viscous  Simulation 


1,000-10,000  time  steps.  A  further  progression  to  large  eddy  simulation  of  complex  configurations  would 
require  even  greater  resources.  The  following  estimate  is  due  to  ^V.H.  Jou  [98].  Suppose  that  a  conservative 
estimate  of  the  size  of  eddies  in  a  boundary  layer  that  ought  to  be  resolved  is  1/5  of  the  boundary  layer 
thickness.  Assuming  that  10  points  are  needed  to  resolve  a  single  eddy,  the  mesh  interval  should  then  be  1/50 
of  the  boundary  layer  thickness.  Moreover,  since  the  eddies  are  three-dimensional,  the  same  mesh  interval 
should  be  Tised  in  all  three  directions.  Now,  if  the  boundary  layer  thickness  is  of  the  order  of  0.01  of  the 
chord  length,  5,000  intervals  will  be  needed  in  the  chordwise  direction,  and  for  a  wing  with  an  aspect  ratio 
of  10,  50,000  intervals  will  be  needed  in  the  spanwise  direction.  Thus,  of  the  order  of  50  x  5, 000  x  50, 000 
or  12.5  billion  mesh  points  would  be  needed  in  the  boundary  layer.  If  the  time  dependent  behavior  of  the 
eddies  is  to  be  fully  resolved  using  time  steps  on  the  order  of  the  time  for  a  wave  to  pass  through  a  mesh 
interval,  and  one  allows  for  a  total  time  equal  to  the  time  required  for  waves  to  travel  three  times  the  length 
of  the  chord,  of  the  order  of  15,000  time  steps  would  be  needed.  A  more  refined  estimate  which  allows  for 
the  varying  thickness  of  the  boundary  layer,  recently  made  by  Spalart  et  al.  [179],  suggests  an  even  more 
severe  requirement.  Performance  beyond  the  terafiop  (10^^  operations  per  second)  will  be  needed  to  attempt 
calculations  of  this  nature,  which  also  have  an  information  content  far  beyond  what  is  needed  for  enginering 
analysis  and  design.  The  designer  does  not  need  to  know  the  details  of  the  eddies  in  the  boundary  layer. 
The  primary  purpose  of  such  calculations  is  to  improve  the  calculation  of  averaged  quantities  such  as  skin 
friction,  and  the  prediction  of  global  behavior  such  as  the  onset  of  separation.  The  main  current  use  of 
Navier-Stokes  and  large  eddy  simulations  is  to  try  to  gain  an  improved  insight  into  the  physics  of  turbulent 
flow,  which  may  in  turn  lead  to  the  development  of  more  comprehensive  and  reliable  turbulence  models. 

2.3.  Turbulence  Modeling.  It  is  doubtful  whether  a  universally  valid  turbulence  model,  capable  of 
describing  all  complex  flows,  could  be  devised  [61].  Algebraic  models  [35,  14]  have  proved  fairly  satisfactory 
for  the  calculation  of  attached  and  slightly  separated  wing  flows.  These  models  rely  on  the  boundary  layer 
concept,  usually  incorporating  separate  formulas  for  the  inner  and  outer  layers,  and  they  require  an  estimate 
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of  a  length  scale  which  depends  on  the  thickness  of  the  boundary  layer.  The  estimation  of  this  quantity  by 
a  search  for  a  maviTmim  of  the  vorticity  times  a  distance  to  the  wall,  as  in  the  Baldwin-Lomax  model,  can 
lead  to  ambiguities  in  internal  flows,  and  also  in  complex  vortical  flows  over  slender  bodies  and  highly  swept 
or  delta  wings  [47,  123].  The  Johnson-King  model  [96],  which  allows  for  non-equilibrium  effects  through  the 
introduction  of  an  ordinary  differential  equation  for  the  maximum  shear  stress,  has  improved  the  prediction 
of  flows  with  shock  induced  separation  [165,  99]. 

Closure  models  depending  on  the  solution  of  transport  equations  are  widely  accepted  for  industrial 
applications.  These  models  eliminate  the  need  to  estimate  a  length  scale  by  detecting  the  edge  of  the 
boundary  layer.  Eddy  viscosity  models  typically  use  two  equations  for  the  turbulent  kinetic  energy  k  and 
the  dissipation  rate  e,  or  a  pair  of  equivalent  quantities  [97, 199, 180, 1, 131, 40] .  Models  of  this  type  generally 
tend  to  present  difficulties  in  the  region  very  close  to  the  wall.  They  also  tend  to  be  badly  conditioned  for 
numerical  solution.  The  k  -  I  model  [173]  is  designed  to  alleviate  this  problem  by  taking  advantage  of  the 
hnear  behaviour  of  the  length  scale  I  near  the  wall.  In  an  alternative  approach  to  the  design  of  models  which 
are  more  amenable  to  munerical  solution,  new  models  requiring  the  solution  of  one  transport  equation  have 
recently  been  introduced  [13,  178].  The  performance  of  the  algebraic  models  remains  competitive  for  wing 
flows,  but  the  one-  and  two-equation  models  show  promise  for  broader  classes  of  flows.  In  order  to  achieve 
greater  iiniversahty,  research  is  also  being  pursued  on  more  complex  Reynolds  stress  transport  models,  which 
require  the  solution  of  a  larger  number  of  transport  equations. 

Another  direction  of  research  is  the  attempt  to  devise  more  rational  models  via  renormalization  group 
(RNG)  theory  [203, 174].  Both  algebraic  and  two-equation  it -e  models  devised  by  this  approach  have  shown 
promising  results  [124], 

The  selection  of  sufficiently  accurate  mathematical  models  and  a  judgment  of  their  cost  effectiveness 
tiltimately  rests  with  industry.  Aircraft  and  spacecraft  designs  normally  pass  through  the  three  phases  of 
conceptual  design,  prehminary  design,  and  detailed  design.  Correspondingly,  the  appropriate  CFD  models 
will  vary  in  complexity.  In  the  conceptual  and  prehnoinary  design  phases,  the  emphasis  wiU  be  on  relatively 
simple  models  which  can  give  results  with  very  rapid  turn-around  and  low  computer  costs,  in  order  to  evaluate 
alternative  configurations  and  perform  quick  parametric  studies.  The  detailed  design  stage  requires  the  most 
complete  simulation  that  can  be  achieved  with  acceptable  cost.  In  the  past,  the  low  level  of  confidence  that 
could  be  placed  on  numerical  predictions  has  forced  the  extensive  use  of  wind  tunnel  testing  at  an  early  stage 
of  the  design.  This  practice  was  very  expensive.  The  hmited  number  of  models  that  could  be  fabricated  also 
hmited  the  range  of  design  variations  that  could  be  evaluated.  It  can  be  anticipated  that  in  the  future,  the 
role  of  wind  tunnel  testing  in  the  design  process  will  be  more  one  of  verification.  Experimental  research  to 
improve  our  understanding  of  the  physics  of  complex  flows  will  continue,  however,  to  play  a  vital  role. 

3.  CFD  Algorithms. 

3.1.  Difficulties  of  Flow  Simulation.  The  computational  simulation  of  fluid  flow  presents  a  number 
of  severe  challenges  for  algorithm  design.  At  the  level  of  inviscid  modeling,  the  inherent  nonlinearity  of  the 
fluid  flow  equations  leads  to  the  formation  of  singularities  such  as  shock  waves  and  contact  discontinuities. 
Moreover,  the  geometric  configurations  of  interest  are  extremely  complex,  and  generally  contain  sharp  edges 
which  lead  to  the  shedding  of  vortex  sheets.  Extreme  gradients  near  stagnation  points  or  wing  tips  may  also 
lead  to  numerical  errors  that  can  have  global  influence.  Numerically  generated  entropy  may  be  convected 
from  the  leading  edge  for  example,  causing  the  formation  of  a  numerically  induced  boundary  layer  which 
can  lead  to  separation.  The  need  to  treat  exterior  domains  of  infinite  extent  is  also  a  source  of  difficulty. 
Boundary  conditions  imposed  at  artificial  outer  boundaries  may  cause  reflected  waves  which  significantly 
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interfere  with  the  flow.  When  viscous  effects  are  also  included  in  the  simulation,  the  extreme  difference  of 
the  scales  in  the  viscous  boundary  layer  and  the  outer  flow,  which  is  essentially  inviscid,  is  another  source 
of  difficulty,  forcing  the  use  of  meshes  with  extreme  variations  in  mesh  interval.  For  these  reasons  CFD,  has 
been  a  driving  force  for  the  development  of  numerical  algorithms. 

3.2.  Structured  and  Unstructured  Meshes.  The  algorithm  designer  faces  a  number  of  critical  deci¬ 
sions.  The  first  choice  that  must  be  made  is  the  nature  of  the  mesh  used  to  divide  the  flow  field  into  discrete 
subdomains.  The  discretization  procedure  must  allow  for  the  treatment  of  complex  configurations.  The 
principal  alternatives  are  Cartesian  meshes,  body-fitted  curvilinear  meshes,  and  unstructured  tetrahedral 
meshes.  Each  of  these  approaches  has  advantages  which  have  led  to  their  use.  The  Cartesian  mesh  mini¬ 
mizes  the  complexity  of  the  algorithm  at  interior  points  and  facilitates  the  use  of  hi^  order  discretization 
procedures,  at  the  expense  of  greater  complexity,  and  possibly  a  loss  of  accuracy,  in  the  treatment  of  bound¬ 
ary  conditions  at  curved  surfaces.  This  difficulty  may  be  alleviated  by  using  mesh  refinement  procedures 
near  the  surface.  With  their  aid,  schemes  which  use  Cartesian  meshes  have  recently  been  developed  to  treat 
very  complex  configurations  [130,  166,  26,  103]. 

Body-fitted  meshes  have  been  widely  used  and  are  particularly  well  suited  to  the  treatment  of  viscous 
flow  because  they  readily  allow  the  mesh  to  be  compressed  near  the  body  surface.  With  this  approach,  the 
problem  of  mesh  generation  itself  has  proved  to  be  a  major  pacing  item.  The  most  commonly  used  procedures 
are  algebraic  transformations  [11,  52,  55, 175],  methods  based  on  the  solution  of  elliptic  equations,  pioneered 
by  Thompson  [189,  190,  176,  177],  and  methods  based  on  the  solution  of  hyperbolic  equations  marching  out 
from  the  body  [181].  In  order  to  treat  very  complex  configurations  it  generally  proves  expedient  to  use  a 
multiblock  [198,  167]  procedure,  with  separately  generated  meshes  in  each  block,  which  may  then  be  patched 
at  block  faces,  or  allowed  to  overlap,  as  in  the  Chimera  scheme  [23,  24].  While  a  number  of  interactive 
software  systems  for  grid  generation  have  been  developed,  such  as  EAGLE,  GRIDGEN,  and  ICEM,  the 
generation  of  a  satisfactory  grid  for  a  very  complex  configuration  may  require  months  of  effort. 

The  alternative  is  to  use  an  unstnictured  mesh  in  which  the  domain  is  subdivided  into  tetrahedra. 
This  in  turn  requires  the  development  of  solution  algorithms  capable  of  yielding  the  required  accuracy  on 
unstructmed  meshes.  This  approach  has  been  gaining  acceptance,  as  it  is  becoming  apparent  that  it  can  lead 
to  a  speed-up  and  reduction  in  the  cost  of  mesh  generation  that  more  than  offsets  the  increased  complexity 
and  cost  of  the  flow  simulations.  Two  competing  procedures  for  generating  triangulations  which  have  both 
proved  successful  are  Delaunay  triangulation  [48,  17],  based  on  concepts  introduced  at  the  beginning  of  the 
century  by  Voronoi  [196],  and  the  moving  front  method  [119]. 


3.3.  Finite  Diflference,  Finite  Volume,  and  Finite  Element  Schemes.  Associated  with  choice 
of  mesh  type  is  the  formulation  of  the  discretization  procedure  for  the  equations  of  fluid  flow,  which  can  be 
expressed  as  differential  conservation  laws.  In  the  Cartesian  tensor  notation,  let  Xj,  be  the  coordinates,  p,  p, 
T,  and  E  the  pressure,  density,  temperature,  and  total  energy,  and  Ui  the  velocity  components.  Using  the 
convention  that  summation  over  j  =  1,2, 3  is  implied  by  a  repeated  subscript  j,  each  conservation  equation 
has  the  form 


(1) 


dw  dfj  df-g.  _  ^ 

dt  dxj  dxj 


where  fj  are  the  inviscid  (convective)  fluxes,  and  fvj  are  the  viscous  fluxes.  For  the  mass  equation 


w  =  p,  fj  =  puj. 
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For  the  i  momentum  equation 


Wi  • —  fij  —  p^i^j  4”  P^ij^  fvij  ^ij'i 


where  aij  is  the  viscous  stress  tensor.  For  the  energy  equation 


w  =  pE,  fj  -  pHuj,  -  ajkUk 


-k 


dT 

dXj’ 


where  k  is  the  coefficient  of  thermal  conductivity.  The  pressure  is  related  to  the  density  and  energy  by  the 
equation  of  state 


(2) 


p=i'y-l)p[E-  -Ui% 


in  which  7  is  the  ratio  of  specific  heats  and  the  stagnation  enthalpy  is  given  by 


H  =  E  +  ?- 
P 


while 


E  —  CxjT  “j” 


where  is  the  specific  heat  at  constant  volume.  In  the  Navier-Stokes  equations  the  viscous  stresses  are 
assumed  to  be  linearly  proportional  to  the  rate  of  strain,  or 


where  ji  and  A  are  the  coefficients  of  viscosity  and  bulk  viscosity,  and  usually  A  =  — 2/i/3.  The  coefficient  of 
thermal  conductivity  and  temperature  are  given  by  the  relations 


k 


rp  _ 

Pr'  Rp 


where  Cp  is  the  specific  heat  at  constant  pressure,  R  is  the  gas  constant,  and  Pr  is  the  Prandtl  number. 

The  finite  difference  method,  which  requires  the  use  of  a  Cartesian  or  a  structured  curvilinear  mesh, 
directly  approximates  the  differential  operators  appearing  in  these  equations.  In  the  finite  volume  method 
[120],  the  discretization  is  accomplished  by  dividing  the  domain  of  the  flow  into  a  large  number  of  small 
subdomains,  and  applying  the  conservation  laws  in  the  integral  form 

4  /  wdV+  f  f-dS  =  0. 

dt  Jn  JdQ 


Here  F  is  the  flux  appearing  in  equation  (1)  and  dS  is  the  directed  surface  element  of  the  boundary  dO,  of 
the  domain  Cl.  The  use  of  the  integral  form  has  the  advantage  that  no  assumption  of  the  differentiability 
of  the  solutions  is  impHed,  with  the  result  that  it  remains  a  valid  statement  for  a  subdomain  containing  a 
shock  wave.  In  general  the  subdomains  could  be  arbitrary,  but  it  is  convenient  to  use  either  hexahedral  cells 
in  a  body  conforming  curvilinear  mesh  or  tetrahedrons  in  an  unstructured  mesh. 

Alternative  discretization  schemes  may  be  obtained  by  storing  flow  variables  at  either  the  cell  centers  or 
the  vertices.  These  variations  are  illustrated  in  Figxire  3  for  the  two-dimensional  case.  With  a  cell-centered 
scheme  the  discrete  conservation  law  takes  the  form 
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Fig.  3.  Sivuctured  and  Unstructured  Discretizations. 


(4) 


^«,V+'£t-S  =  0, 

faces 


where  V  is  the  cell  volume,  and  f  is  now  a  numerical  estimate  of  the  flux  vector  through  each  face,  f  may 
be  evaluated  from  values  of  the  flow  variables  in  the  cells  separated  by  each  face,  using  upwind  biasing  to 
allow  for  the  directions  of  wave  propagation.  With  hexahedral  cells,  equation  (4)  is  very  similar  to  a  finite 
difference  scheme  in  curvilinear  coordinates.  Under  a  transformation  to  curvilinear  coordinates  ,  equation 
(1)  becomes 


(5)  ■ 


dt 


+ w< 


=  0, 


where  J  is  the  Jacobian  determinant  of  the  transformation  matrix  .  The  transformed  flux  J 
corresponds  to  the  dot  product  of  the  flux  f  with  a  vector  face  area  while  J  represents  the  transfor¬ 

mation  of  the  cell  volume.  The  finite  volume  form  (4)  has  the  advantages  that  it  is  valid  for  both  structured 
and  unstructured  meshes,  and  that  it  assures  that  a  uniform  flow  exactly  satisfies  the  equations,  because 
S  =  0  for  a  closed  control  volume.  Finite  difference  schemes  do  not  necessarily  satisfy  this  constraint 
because  of  the  discretization  errors  in  evaluating  ^  and  the  inversion  of  the  transformation  matrix.  A 
cell- vertex  finite  volume  scheme  can  be  derived  by  taking  the  union  of  the  cells  surrounding  a  given  vertex 
as  the  control  volume  for  that  vertex  [64,  81,  156].  In  equation  (4),  V  is  now  the  sum  of  the  volumes  of  the 
surrounding  cells,  while  the  flux  balance  is  evaluated  over  the  outer  faces  of  the  polyhedral  control  volume. 
In  the  absence  of  upwind  biasing  the  flux  vector  is  evaluated  by  averaging  over  the  corners  of  each  face.  This 
has  the  advantage  of  remaining  accurate  on  an  irregular  or  unstructured  mesh.  An  alternative  route  to  the 
discrete  equations  is  provided  by  the  finite  element  method.  Whereas  the  finite  difference  and  finite  volume 
methods  approximate  the  differential  and  integral  operators,  the  finite  element  method  proceeds  by  inserting 
an  approximate  solution  into  the  exact  equations.  On  multiplying  by  a  test  fimction  (j>  and  integrating  by 
parts  over  space,  one  obtains  the  weak  form 
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which  is  also  valid  in  the  presence  of  discontinuities  in  the  flow.  In  the  Galerkin  method  the  approximate 
solution  is  expanded  in  terms  of  the  same  family  of  functions  as  those  from  which  the  test  functions  are 
drawn.  By  choosing  test  functions  with  local  support,  separate  equations  are  obtained  for  each  node.  For 
example,  if  a  tetrahedral  mesh  is  used,  and  (j>  is  piecewise  linear,  with  a  nonzero  value  only  at  a  single  node, 
the  equations  at  each  node  have  a  stencil  which  contains  only  the  nearest  neighbors.  In  this  case  the  finite 
element  approximation  corresponds  closely  to  a  finite  volume  scheme.  If  a  piecewise  linear  approximation 
to  the  flux  f  is  used  in  the  evaluation  of  the  integrals  on  the  right  hand  side  of  equation  (6),  these  integrals 
reduce  to  formulas  which  are  identical  to  the  flux  balance  of  the  finite  volume  scheme. 

Thus  the  finite  difference  and  finite  volume  methods  lead  to  essentially  similar  schemes  on  structured 
meshes,  while  the  finite  volume  method  is  essentially  equivalent  to  a  finite  element  method  with  linear 
elements  when  a  tetrahedral  mesh  is  used.  Provided  that  the  flow  equations  are  expressed  in  the  conservation 
law  form  (1),  all  three  methods  lead  to  an  exact  cancellation  of  the  fluxes  through  interior  cell  boundaries,  so 
that  the  conservative  property  of  the  equations  is  preserved.  The  important  role  of  this  property  in  ensuring 
correct  shock  jump  conditions  was  pointed  out  by  Lax  and  Wendroff  [105]. 


3.4.  Non-oscillatory  Shock  Capturing  Schemes. 

3.4.1.  Local  Extremum  Diminishing  (LED)  Schemes.  The  discretization  procedures  which  have 
been  described  in  the  last  section  lead  to  nondissipative  approximations  to  the  Euler  equations.  Dissipative 
terms  may  be  needed  for  two  reasons.  The  first  is  the  possibility  of  undamped  oscillatory  modes.  The  second 
reason  is  the  need  for  the  clean  capture  of  shock  waves  and  contact  discontinuities  without  undesirable 
oscillations.  An  extreme  overshoot  could  result  in  a  negative  value  of  an  inherently  positive  quantity  such  as 
the  pressure  or  density.  The  next  sections  summarize  a  unified  approach  to  the  construction  of  nonoscillatory 
schemes  via  the  introduction  of  controlled  diffusive  and  antidiffusive  terms.  This  is  the  line  adhered  to  in 
the  author’s  own  work. 

The  development  of  non-oscillatory  schemes  has  been  a  prime  focus  of  algorithm  research  for  compressible 
flow.  Consider  a  general  semi-discrete  scheme  of  the  form 


(7) 


-Vj). 


A  maximum  cannot  increase  and  a  minimum  cannot  decrease  if  the  coefficients  Cjk  are  non-negative,  since 
at  a  maximum  <  0,  and  at  a  minimum  Vk  ^  Vj  '>  0,  Thus  the  condition 

(8)  Cjk  >0,  k^j 


is  sufficient  to  ensure  stability  in  the  maximum  norm.  Moreover,  if  the  scheme  has  a  compact  stencil,  so  that 
Cjk  =  0  when  j  and  k  are  not  nearest  neighbors,  a  local  maximum  cannot  increase  and  local  minimum  cannot 
decrease.  This  local  extremum  diminishing  (LED)  property  prevents  the  birth  and  growth  of  oscillations. 
The  one-dimensional  conservation  law 


du 

m 


+  =  0 


provides  a  useful  model  for  analysis.  In  this  case  waves  are  propagated  with  a  speed  a{u)  =  |£,  and  the 
solution  is  constant  along  the  characteristics  ^  —  q>{u).  Thus  the  LED  property  is  satisfied.  In  fact  the 
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total  variation 


dx 


of  a  solution  of  this  equation  does  not  increase,  provided  that  any  discontinuity  appearing  in  the  solution 
satisfies  an  entropy  condition  [104].  Harten  proposed  that  difference  schemes  ought  to  be  designed  so  that  the 
discrete  total  variation  cannot  increase  [65],  If  the  end  values  are  fixed,  the  total  variation  can  be  expressed 
as 


TV{u)  =  2i^ 


maxima 


E 


minima 


)■ 


Thus  a  LED  scheme  is  also  total  variation  diminishing  (TVD).  Positivity  conditions  of  the  type  expressed  in 
equations  (7)  and  (8)  lead  to  diagonally  dominant  schemes,  and  are  the  key  to  the  elimination  of  improper 
oscillations.  The  positivity  conditions  may  be  realized  by  the  introduction  of  diffusive  terms  or  by  the  use  of 
upwind  biasing  in  the  discrete  scheme.  Unfortunately,  they  may  also  lead  to  severe  restrictions  on  accuracy 
imless  the  coefficients  have  a  complex  nonlinear  dependence  on  the  solution. 


3.4.2.  Artificial  Diffusion  and  Upwinding.  Following  the  pioneering  work  of  Godunov  [60] ,  a  variety 
of  dissipative  and  upwind  schemes  designed  to  have  good  shock  capturing  properties  have  been  developed 
during  the  past  two  decades  [182,  27,  107,  108,  163,  142,  65,  141,  185,  8,  79,  204,  73,  201,  16,  15,  17].  If  the 
one-dimensional  scalar  conservation  law 


(9) 


dv  d  .  . 
dt  ' 


is  represented  by  a  three  point  scheme 

^  («i+i  -  ~  ’ 


the  scheme  is  LED  if 


(10) 


c+  ,  >  0,  c7  ,  >  0. 
3+h  -  3-h 


A  conservative  semidiscrete  approximation  to  the  one-dimensional  conservation  law  can  be  derived  by  sub¬ 
dividing  the  line  into  cells.  Then  the  evolution  of  the  value  vj  in  the  jth  ceU  is  given  by 


(11) 


A  dvj 


hA_i 
3  2 


0, 


where  h,.i  is  an  estimate  of  the  flux  between  cells  j  and  j  +  1.  The  simplest  estimate  is  the  arithmetic 

J-r  2 

average  {fj+i  +  fj)/2,  but  this  leads  to  a  scheme  that  does  not  satisfy  the  positivity  conditions.  To  correct 
this,  one  may  add  a  dissipative  term  and  set 

(12)  hj+i  =  i  (/j+i  +  fj)  -  1  (uj+i  -  Vj) . 


In  order  to  estimate  the  required  value  of  the  coefficient  Q!j+ 1 ,  let  aj+i  be  a  numerical  estimate  of  the  wave 
speed 


(13) 


^J+i  = 


§77^  if  Vj+i^Vj 
%  if  Vj+i  =  Vj 

V=Vj 
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Then 


where 


^'^j+h~'^j+'>-  ‘^i> 


and  the  LED  condition  (10)  is  satisfied  if 
(14) 


^  2 


“i+l 


If  one  takes 


S+2 


one  obtains  the  first  order  upwind  scheme 


f  fj  if  flj+i  >  0 
1  fj+i  if  “j+i  <  0 


This  is  the  least  diffusive  first  order  scheme  which  satisfies  the  LED  condition.  In  this  sense  upwinding 
is  a  natural  approach  to  the  construction  of  non-oscillatory  schemes.  It  may  be  noted  that  the  successful 
treatment  of  transonic  potential  flow  also  involved  the  use  of  upwind  biasing.  This  was  first  introduced  by 
Murman  and  Cole  to  treat  the  transonic  small  disturbance  equation  [135]. 

Another  important  requirement  of  discrete  schemes  is  that  they  should  exclude  nonphysical  solutions 
which  do  not  satisfy  appropriate  entropy  conditions  [106],  which  require  the  convergence  of  characteristics 
towards  admissible  discontinuities.  This  places  more  stringent  boimds  on  the  minimum  level  of  numerical 
viscosity  [121, 188,  140,  143].  In  the  case  that  the  numerical  flux  function  is  strictly  convex.  Also  has  recently 
proved  [2]  that  it  is  sufficient  that 


>  max 


,esign(Uj+i -Uj)| 


for  e  >  0.  Thus  the  numerical  viscosity  should  be  rounded  out  and  not  allowed  to  reach  zero  at  a  point 
where  the  wave  speed  a(u)  =  approaches  zero.  This  justifies,  for  example,  Harten’s  entropy  fix  [65]. 

Higher  order  schemes  can  be  constructed  by  introducing  higher  order  diffusive  terms.  Unfortunately 
these  have  larger  stencils  and  coefficients  of  varying  sign  which  are  not  compatible  with  the  conditions  (8) 
for  a  LED  scheme,  and  it  is  known  that  schemes  which  satisfy  these  conditions  are  at  best  first  order  accurate 
in  the  neighborhood  of  an  extremum.  It  proves  useful  in  the  following  development  to  introduce  the  concept 
of  essentially  local  extremum  diminishing  (ELED)  schemes.  These  are  defined  to  be  schemes  which  satisfy 
the  condition  that  in  the  limit  as  the  mesh  width  Ax  —>■  0,  local  maxima  are  non-increasing,  and  local 
minima  are  non-decreasing. 

3.4.3.  High  Resolution  Switched  Schemes:  Jameson-Schmidt-Turkel  (JST)  Scheme.  Higher 
order  non-oscillatory  schemes  can  be  derived  by  introducing  anti-diffusive  terms  in  a  controlled  manner. 
An  early  attempt  to  produce  a  high  resolution  scheme  by  this  approach  is  the  Jameson-Schmidt-Turkel 
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(JST)  scheme  [93].  Suppose  that  anti-diffusive  terms  are  introduced  by  subtracting  neighboring  differences 
to  produce  a  third  order  diffusive  flux 

(15)  I  }  - 

which  is  an  approximation  to  ^aAx^^.  The  positivity  condition  (8)  is  violated  by  this  scheme.  It  proves 
that  it  generates  substantial  oscillations  in  the  vicinity  of  shock  waves,  which  can  be  eliminated  by  switching 
locally  to  the  first  order  scheme.  The  JST  scheme  therefore  introduces  blended  diffusion  of  the  form 


(16) 


2 

J  + 


i  (' 


2Au^+i  + 


The  idea  is  to  use  variable  coefficients  and  which  produce  a  low  level  of  diffusion  in  regions  where 
the  solution  is  smooth,  but  prevent  oscillations  near  discontinuities.  If  e^.^i  is  constructed  so  that  it  is  of 
order  where  the  solution  is  smooth,  while  is  of  order  unity,  both  terms  in  dj+i  will  be  of  order 

The  JST  scheme  has  proved  very  effective  in  practice  in  numerous  calculations  of  complex  steady  flows, 
and  conditions  under  which  it  could  be  a  total  variation  diminishing  (TVD)  scheme  have  been  examined 

/o'! 

by  Swanson  and  Turkel  [184].  An  alternative  statement  of  sufficient  conditions  on  the  coefficients  and 
for  the  JST  scheme  to  be  LED  is  as  follows: 

2 

Theorem  3,1  (Positivity  of  the  JST  scheme).  Suppose  that  whenever  either  orvj  is  an  extremum 
the  coefficients  of  the  JST  scheme  satisfy 


(17) 


=  0. 


Then  the  JST  scheme  is  local  extremum  diminishing  (LED). 

Proof:  We  need  only  consider  the  rate  of  change  ofv  at  extremal  points.  Suppose  thatvj  is  an  extremum. 
Then 


M)  _,(4)  _o 


and  the  semi-discrete  scheme  (11)  reduces  to 

A 


and  each  coefficient  has  the  required  sign.  □ 

.(2)  .(4) 


J2) 


ax.L  Av. 


( Avj_i, 

\3-2  2  2  y  J  2 


In  order  to  construct  ev^  i  and  e\  ^  i  with  the  desired  properties  define 
J”2  3  2 


(18) 


R{u,v)  = 


r  1 

u—v 

I 

1 

H+hl  1 

1 

iiu^O  OT  V  ^0 
0  if  u  =  i;  =  0, 


where  g  is  a  positive  integer.  Then  R{u,v)  =  lifu  and  v  have  opposite  signs.  Otherwise  R{u,v)  <  1.  Now 
set 


Qj  =  1  =  max(Qj, Q,-+i). 
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and 

(19) 


=  i. 


3.4.4,  Symmetric  Limited  Positive  (SLIP)  Scheme.  An  alternative  route  to  high  resolution  with¬ 
out  oscillation  is  to  introduce  flux  limiters  to  guarantee  the  satisfaction  of  the  positivity  condition  (8).  The 
use  of  limiters  dates  back  to  the  work  of  Boris  and  Book  [27],  A  particularly  simple  way  to  introduce  lim¬ 
iters,  proposed  by  the  author  in  1984  [79],  is  to  use  flux  limited  dissipation.  In  this  scheme  the  third  order 
diffusion  defined  by  equation  (15)  is  modified  by  the  insertion  of  limiters  which  produce  an  equivalent  three 
point  scheme  with  positive  coefficients.  The  original  scheme  [79]  can  be  improved  in  the  following  manner  so 
that  less  restrictive  flux  limiters  are  required.  Let  L(u,u)  be  a  limited  average  of  u  and  v  with  the  following 
properties: 

PL  L{%v)  =  L{v,u) 

P2.  L{au^av)  =  aL{u,v) 

P3.  L{u,u)  =  u 

P4.  L{ujV)  =  0  if  n  and  v  have  opposite  signs:  otherwise  L(u,  u)  has  the  same  sign  as  u  and  v. 
Properties  (P1-P3)  are  natural  properties  of  an  average.  Property  (P4)  is  needed  for  the  construction  of  a 
LED  or  TVD  scheme. 

It  is  convenient  to  introduce  the  notation 

<Pir)  =  L{l,r)  =  L{r,l), 

where  according  to  (P4)  (f>{r)  >  0.  It  follows  from  (P2)  on  setting  a  =  ^  or  ^  that 

Also  it  follows  on  setting  v  —  1  and  u  =  r  that 

4>{r)  =  r4>  ^1)  . 

Thus,  if  there  exists  r  <  0  for  which  ^(r)  >  0,  then  (p  (i)  <  0.  The  only  way  to  ensure  that  ^(r)  >  0  is  to 
require  ^(r)  =  0  for  all  r  <  0,  corresponding  to  property  (P4). 

Now  one  defines  the  diffusive  flux  for  a  scalar  conservation  law  as 


(20) 

Set 

+  _ 

1 

<1 

1 

V  — 

and 

L{Avj+^,Avj_i)  =  <t>{r+)Avj_x 
L{Avj_3,Avj+i)  =  ^(r-)Avj+L. 

Then, 

(21)  i  <^(r+)|  Av^_i. 
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Thus  the  scheme  satisfies  the  LED  condition  if  ^  for  all  j,  and  (I>{t)  >  0,  which  is  assured 

by  property  (P4)  on  L.  At  the  same  time  it  follows  from  property  (P3)  that  the  first  order  diffusive  flux  is 
canceled  when  Av  is  smoothly  varying  and  of  constant  sign.  Schemes  constructed  by  this  formulation  will 
be  referred  to  as  symmetric  limited  positive  (SLIP)  schemes.  This  result  may  be  summarized  as 

Theorem  3.2  (Positivity  of  the  SLIP  scheme).  Suppose  that  the  discrete  conservation  law  (11)  contains 
a  limited  diffusive  flux  as  defined  by  equation  (^0).  Then  the  positivity  condition  (14)}  together  with  the 
properties  (P1~P4)  for  limited  averages,  are  sufficient  to  ensure  satisfaction  of  the  LED  principle  that  a  local 
maximum  cannot  increase  and  a  local  minimum  cannot  decrease.  □ 

A  variety  of  limiters  may  be  defined  which  meet  the  requirements  of  properties  (P1-P4).  Define 

S{u,v)  =  i  {sign(u)  +  sign(u)} 

which  vanishes  is  u  and  v  have  opposite  signs. 

Then  two  limiters  which  are  appropriate  are  the  following  well-known  schemes: 

1.  Minmod: 

L{u,v)  =  S'(u,u)min(|u|,  lu|) 

2.  Van  Leer: 

In  order  to  produce  a  family  of  limiters  which  contains  these  as  special  cases  it  is  convenient  to  set 

L{u,v)  =  ^D{u,v){u  +  v), 

where  D{u,v)  is  a  factor  which  should  deflate  the  arithmetic  average,  and  become  zero  if  v,  and  v  have 
opposite  signs.  Take 

(22)  D{u,v)  =  1-  R{u,v)  =  1  - 

where  i?(u,  v)  is  the  same  function  that  was  introduced  in  the  JST  scheme,  and  g  is  a  positive  integer.  Then 
D{u,v)  =  0  if  ix  and  u  have  opposite  signs.  Also  if  g  =  1,  L{u,v)  reduces  to  minmod,  while  if  g  =  2,  L{u,v) 
is  equivalent  to  Van  Leer’s  limiter.  By  increasing  q  one  can  generate  a  sequence  of  limited  averages  which 
approach  a  limit  defined  by  the  arithmetic  mean  truncated  to  zero  when  u  and  v  have  opposite  signs. 

When  the  terms  are  regrouped,  it  can  be  seen  that  with  this  limiter  the  SLIP  scheme  is  exactly  equivalent 
to  the  JST  scheme,  with  the  switch  is  defined  as 


This  formulation  thus  unifies  the  JST  and  SLIP  schemes. 

3.4.5.  Essentially  Local  Extremum  Diminishing  (ELED)  Scheme  with  Soft  Limiter.  The 

limiters  defined  by  the  formula  (22)  have  the  disadvantage  that  they  are  active  at  a  smooth  extrema,  reducing 
the  local  accuracy  of  the  scheme  to  first  order.  In  order  to  prevent  this,  the  SLIP  scheme  can  be  relaxed 


u  —  v 
|u|  +  \v 
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to  give  an  essentially  local  extremnm  diminishing  (ELED)  scheme  which  is  second  order  accurate  at  smooth 
extrema  by  the  introduction  of  a  threshold  in  the  limited  average.  Therefore  redefine  D{u,  v)  as 

q 

r'C'UjV)  1  maxdul  +  |u| ,  eAa;’’)  ’ 

where  r-  =  |,  g  >  2.  This  reduces  to  the  previous  definition  if  |u|  +  |u|  >  eAx^ . 

In  any  region  where  the  solution  is  smooth,  —  Avj_i  is  of  order  Ax^ .  In  fact  if  there  is  a  smooth 

extremum  in  the  neighborhood  of  Vj  or  Uj+i,  a  Taylor  series  expansion  indicates  that  Avj^s,  Avj^^  and 
AVj_i  are  each  individually  of  order  Ax^,  since  ^  =  0  at  the  extremum.  It  may  be  verified  that  second 
order  accuracy  is  preserved  at  a  smooth  extremum  if  g  >  2.  On  the  other  hand  the  limiter  acts  in  the  usual 
way  if  j  A'tk  I  3  or  AVj^  s  ^  cAx^ ,  and  it  may  also  be  verified  that  in  the  limit  Ax  ^  0  local  maxima 
are  non  increasing  and  local  minima  are  non  decreasing  [86].  Thus  the  scheme  is  essentially  local  extremum 
diminishing  (ELED). 

The  effect  of  the  “soft  limiter”  is  not  only  to  improve  the  accuracy:  the  introduction  of  a  threshold 
below  which  extrema  of  small  amplitude  are  accepted  also  usually  results  in  a  faster  rate  of  convergence  to 
a  steady  state,  and  decreases  the  likelyhood  of  limit  cycles  in  which  the  limiter  interacts  unfavorably  with 
the  corrections  produced  by  the  updating  scheme.  In  a  scheme  recently  proposed  by  Venkatakrishnan  a 
threshold  is  introduced  precisely  for  this  purpose  [193]. 


3.4.6.  Upstream  Limited  Positive  (USLIP)  Schemes.  By  adding  the  anti-diffusive  correction 
purely  from  the  upstream  side  one  may  derive  a  family  of  upstream  limited  positive  (USLIP)  schemes. 
Corresponding  to  the  original  SLIP  scheme  defined  by  equation  (20),  a  USLIP  scheme  is  obtained  by  setting 

=  “i+i  ^  } 

if  0^+^  >  0,  or 


if  Oj+i  <  0.  If 

form.  Consider  the  case  that  dj^i 


me  recovers  a  standard  high  resolution  upwind  scheme  in  semi-discrete 


>  0  and  a.  1  >  0.  If  one  sets 

J  2 


Av 


Av 


Av, 


r  = 


J  2 


Av,- 


-I 


the  scheme  reduces  to 

i  {<^(r+)a,+.l  +  (2  }  Avj_x. 

To  assure  the  correct  sign  to  satisfy  the  LED  criterion  the  flux  limiter  must  now  satisfy  the  additional 
constraint  that  <^(r)  <  2. 

The  USLIP  formulation  is  essentially  equivalent  to  standard  upwind  schemes  [142,  185].  Both  the  SLIP 
and  USLIP  constructions  can  be  implemented  on  unstructured  meshes  [85,  86].  The  anti-diffusive  terms  are 
then  calculated  by  taking  the  scalar  product  of  the  vectors  defining  an  edge  with  the  gradient  in  the  adjacent 
upstream  and  downstream  cells. 


3.4.7.  Systems  of  Conservation  Laws:  Flux  Splitting  and  Flux-Difference  Splitting.  Steger 
and  Warming  [182]  first  showed  how  to  generalize  the  concept  of  upwinding  to  the  system  of  conservation 
laws 


(24) 


dw  d  .  , 

_  +  _/(„)  =  0 


15 


by  the  concept  of  flux  splitting.  Suppose  that  the  flux  is  split  as  /  =  /+  +  /  where  ^  and  ^  have 
positive  and  negative  eigenvalues.  Then  the  first  order  upwind  scheme  is  produced  by  taking  the  numerical 
flux  to  be 

*6+1* 

This  can  be  expressed  in  viscosity  form  as 

ifj+i  +  //)  +  2  ) 

=  2  (/j+i  +  /?■)“ 

where  the  diffusive  flux  is 

(25)  dj+i  =  2^(/  )i+5' 

Roe  derived  the  alternative  formulation  of  flux  difference  splitting  [163]  by  distributing  the  corrections  due 
to  the  flux  difference  in  each  interval  upwind  and  downwind  to  obtain 

+  (/,■+!  -  /,)-  +  {fj  -  /,-i)+  =  0, 

where  now  the  flux  difference  fj+i  -  fj  is  split.  The  corresponding  diffusive  flux  is 

Following  Roe’s  derivation,  let  Aj^i  be  a  mean  value  Jacobian  matrix  exactly  satisfying  the  condition 

(26)  fj+i  -  fj  =  Aj+i{wj+i  -  Wj). 

Aj^i  may  be  calculated  by  substituting  the  weighted  averages 

■y/ft+Twi+i  +  \/pI^j  „  _  +  y/pjHi 

into  the  standard  formulas  for  the  Jacobian  matrix  A  =  ^.  A  splitting  according  to  characteristic  fields  is 
now  obtained  by  decomposing  as 

(28)  Aj+^=TAT-\ 

where  the  columns  of  T  are  the  eigenvectors  of  Aj^\,  and  A  is  a  diagonal  matrix  of  the  eigenvalues.  Now 
the  corresponding  diffusive  flux  is 

^|Aj+i|(Wj+1  -Wj), 

where 

A,.+  i|=T|A|T-i 

and  |A|  is  the  diagonal  matrix  containing  the  absohite  values  of  the  eigenvalues. 
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3.4.8.  Alternative  Splittings.  Characteristic  splitting  has  the  advantages  that  it  introduces  the  min¬ 
imum  amount  of  diffusion  to  exclude  the  growth  of  local  extrema  of  the  characteristic  variables,  and  that 
with  the  Roe  linearization  it  allows  a  discrete  shock  structure  with  s  single  interior  point.  To  reduce  the 
computational  complexity  one  may  replace  |.4|  by  ocl  where  if  a  is  at  least  equal  to  the  spectral  radius 
maxlA(.4)|,  then  the  positivity  conditions  will  still  be  satisfied.  Then  the  first  order  scheme  simply  has  the 
scalar  diffusive  flux 


(29) 


di4-i  = 


The  JST  scheme  with  scalar  diffusive  flux  captures  shock  waves  with  about  3  interior  points,  and  it  has  been 
widely  used  for  transonic  flow  calculations  because  it  is  both  robust  and  computationally  inexpensive. 

An  intermediate  class  of  schemes  can  be  formulated  by  defining  the  first  order  diffusive  flux  as  a  combi¬ 
nation  of  differences  of  the  state  and  flux  vectors 


(30)  =  ia*+ic(wj+i  Wj)  +  i  (/j+i  fj) , 

where  the  factor  c  is  included  in  the  first  term  to  make  and  dimensionless.  Schemes  of  this  class 
are  fully  upwind  in  supersonic  flow  if  one  takes  —  0  and  ^j+i  —  sign(M)  when  the  absolute  value  of 
the  Mach  number  M  exceeds  1.  The  flux  vector  /  can  be  decomposed  as 


(31) 


f  =  uw  +  fp, 


where 


(32) 


fp  — 


(:v 


V  / 


Then 

(33)  /i+i  fj  =  u(wj+i  wj)+w(uj+i  %■)  +  fp  j+i  fpj’ 
where  u  and  u>  are  the  arithmetic  averages 

«  =  5  (“i+i  +  %■)  >  ^  • 

Thus  these  schemes  are  closely  related  to  schemes  which  introduce  separate  spfittings  of  the  convective  and 
pressure  terms,  sxich  as  the  wave-particle  scheme  [158,  12],  the  advection  upwind  sphtting  method  (AUSM) 
[114,  197],  and  the  convective  upwind  and  split  pressure  (CUSP)  schemes  [85]. 

In  order  to  examine  the  shock  capturing  properties  of  these  various  schemes,  consider  the  general  case 
of  a  first  order  diffusive  flux  of  the  form 

(34)  dj+ 1  =  1  {Wj+i  -wj), 

where  the  matrix  B..i  determines  the  properties  of  the  scheme  and  the  scaling  factor  a,.i  is  included  for 
convenience.  All  the  previous  schemes  can  be  obtained  by  representing  as  a  polynomial  in  the  matrix 

defined  by  equation  (26).  Schemes  of  this  class  were  considered  by  Van  Leer  [107].  According  to  the 
Cayley-Hamilton  theorem,  a  matrix  satisfies  its  own  characteristic  equation.  Therefore  the  third  and  higher 
powers  of  A  can  be  ehminated,  and  there  is  no  loss  of  generality  in  limiting  Bj+i  to  a  polynomial  of  degree 
2, 

(35)  Bj+L  =  aol  +  +  “2^1^  i- 
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The  characteristic  upwind  scheme  for  which  -Bj+i  =  | 

\  =  TJS?T~'^.  Then  aoi  «i)  and  a2  are  determined 

2 


1  is  obtained  by  substituting  =  TAT 

from  the  three  equations 


1 


0^0  +  OiiXk  +  —  \^k  I  J  --  1, 2, 3. 


The  same  representation  remains  valid  for  three  dimensional  flow  because  Aj_^  i  still  has  only  three  distinct 
eigenvalues  u,u  +  c,u  —  c. 


Fig.  4.  Shock  structure  for  single  interior  point. 

3.4.9.  Analysis  of  Stationary  Discrete  Shocks.  The  ideal  model  of  a  discrete  shock  is  illustrated 
in  figure  (4).  Suppose  that  wl  and  wr  are  left  and  right  states  which  satisfy  the  jump  conditions  for  a 
stationary  shock,  and  that  the  corresponding  fluxes  are  fi  =  f{wi.)  and  fn  =  /{wr).  Since  the  shodc  is 
stationary  /l  =  fn-  The  ideal  discrete  shock  has  constant  states  wl  to  the  left  and  wr  to  the  right,  and  a 
single  point  with  an  intermediate  value  wj^.  The  intermediate  value  is  needed  to  allow  the  discrete  solution 
to  correspond  to  a  true  solution  in  which  the  shock  wave  does  not  coincide  with  an  interface  between  two 
mesh  cells. 

Schemes  corresponding  to  one,  two  or  three  terms  in  equation  (35)  are  examined  in  [87].  The  analysis  of 
these  three  cases  shows  that  a  discrete  shock  structure  with  a  single  interior  point  is  supported  by  artificial 
diffusion  that  satisfies  the  two  conditions  that 

1.  it  produces  an  upwind  flux  if  the  flow  is  determined  to  be  supersonic  through  the  interface 

2.  it  satisfies  a  generalized  eigenvalue  problem  for  the  exit  from  the  shock  of  the  form 

(36)  {Aar  -  ocarBar)  {wr  -  Wyi)  =  0, 

where  Aar  is  the  flnearized  Jacobian  matrix  and  Bar  is  the  matrix  defining  the  diffusion  for  the  interface 
AR.  This  follows  from  the  equihbrium  condition  Hra  =  hRR  for  the  cell  j  +  1  in  figure  4.  These  two 
conditions  are  satisfied  by  both  the  characteristic  scheme  and  also  the  CUSP  scheme,  provided  that  the 
coefficients  of  convective  diffusion  and  pressure  differences  are  correctly  balanced.  Scalar  diffusion  does  not 
satisfy  the  first  condition.  In  the  case  of  the  CUSP  scheme  (30)  equation  (36)  reduces  to 

(ara  +  -wa)  =  0 

Thus  WR  -  WA  is  an  eigenvector  of  the  Roe  matrix  Ara,  and  is  the  corresponding  eigenvalue.  Since 

the  eigenvalues  are  u,  u  +  c,  and  u  —  c,  the  only  choice  which  leads  to  positive  diffusion  when  u  >  0  is  u  —  c, 
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yielding  the  relationship 


a*c  =  (1  +  /?)(c- ti),0  <u<c 

This  leads  to  a  one  parameter  family  of  schemes  which  support  the  ideal  shock  structure.  The  term  /3(/i?,-/A) 
contributes  to  the  diffusion  of  the  convective  terms.  Allowing  for  the  split  (31),  the  total  effective  coefficient 
of  convective  diffusion  is  ac  =  o*c  +  /3u.  A  CUSP  scheme  with  low  numerical  diffusion  is  then  obtained  by 
taking  a  =  \M\,  leading  to  the  coefficients  illustrated  in  figure  5 


Fig.  6.  Diffusion  Coefficients. 


3.4.10.  CUSP  and  Characteristic  Schemes  Admitting  Constant  Total  Enthalpy  in  Steady 
Plow.  In  steady  flow  the  stagnation  enthalpy  H  is  constant,  corresponding  to  the  fact  that  the  energy  and 
mass  conservation  equations  are  consistent  when  the  constant  factor  H  is  removed  from  the  energy  equation. 
Discrete  and  semi-discrete  schemes  do  not  necessarily  satisfy  this  property.  In  the  case  of  a  semi-discrete 
scheme  expressed  in  viscosity  form,  equations  (11)  and  (12),  a  solution  with  constant  H  is  admitted  if  the 
viscosity  for  the  energy  equation  reduces  to  the  viscosity  for  the  continuity  equation  with  p  replaced  by 
pH.  When  the  standard  characteristic  decomposition  (28)  is  used,  the  viscous  fluxes  for  p  and  pH  which 
result  from  composition  of  the  fluxes  for  the  characteristic  variables  do  not  have  this  property,  and  H  is  not 
constant  in  the  discrete  solution.  In  practice  there  is  an  excursion  of  H  in  the  discrete  shock  structure  which 
represents  a  local  heat  source.  In  very  high  speed  flows  the  corresponding  error  in  the  temperature  may  lead 
to  a  wrong  prediction  of  associated  effects  such  as  chemical  reactions. 

The  source  of  the  error  in  the  stagnation  enthalpy  is  the  discrepancy  between  the  convective  terms 


u 


P  \ 
pu 

\PH  J 


in  the  flux  vector,  which  contain  pH,  and  the  state  vector  which  contains  pE.  This  may  be  remedied  by 
introducing  a  modified  state  vector 


Then  one  introduces  the  linearization 


/r  -  fL  =  -AhiWhR  -  Whi, ). 


Here  Ah  may  be  calculated  in  the  same  way  as  the  standard  Roe  linearization.  Introduce  the  weighted 
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averages  defined  by  equation  (27).  Then 


/  0 


Ah  = 


I 

7  2 
\  -uH 


1  0  \ 

"r»  ^ 

H  u  J 


The  eigenvalues  of  Ah  are  u,  A"*"  and  A  where 


(37) 


A±  =  ^u± 

27 


— 
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Now  both  CUSP  and  characteristic  schemes  which  preserve  constant  stagnation  enthalpy  in  steady  flow  can 
be  constructed  from  the  modified  Jacobian  matrix  Ah  [87].  These  schemes  also  produce  a  discrete  shock 
structure  with  one  interior  point  in  steady  flow.  Then  one  arrives  at  four  variations  with  this  property,  which 
can  conveniently  be  distinguished  as  the  E-  and  H-CUSP  schemes,  and  the  E-  and  H-characteristic  schemes. 

3.4.11.  High  Order  Godunov  Schemes,  and  Kinetic  Flux  Splitting.  Some  of  the  most  impres¬ 
sive  simulations  of  time  dependent  flows  with  strong  shock  waves  have  been  achieved  with  higher  order 
Godunov  schemes  [201].  In  these  schemes  the  average  value  in  each  cell  is  updated  by  applying  the  inte¬ 
gral  conservation  law  using  interface  fluxes  predicted  from  the  exact  or  approximate  solution  of  a  Riemann 
problem  between  adjacent  cells.  A  higher  order  estimate  of  the  solution  is  then  reconstructed  from  the 
cell  averages,  and  slope  limiters  are  apphed  to  the  reconstruction.  An  example  is  the  class  of  essentially 
non-oscillatory  (ENO)  sdiemes,  which  can  attain  a  very  high  order  of  accuracy  at  the  cost  of  a  substan¬ 
tial  increase  in  computational  complexity  [37,  170,  168,  169].  Methods  based  on  reconstruction  can  also 
be  implemented  on  unstnictured  meshes  [16,  15].  Recently  there  has  been  an  increasing  interest  in  kinetic 
flux  splitting  schemes,  which  use  solutions  of  the  Boltzmann  equation  or  the  BGK  equation  to  predict  the 
interface  fluxes  [50,  41,  54,  153,  202]. 

3.5.  Multidimensional  Schemes.  The  simplest  approach  to  the  treatment  of  multi-dimensional  prob¬ 
lems  on  structured  meshes  is  to  apply  the  one-dimensional  construction  separately  in  each  mesh  direction.  On 
triangulated  meshes  in  two  or  three  dimensions  the  SLIP  and  USLIP  constructions  may  also  be  implemented 
along  the  mf>.<!Vi  edges  [86].  A  substantial  body  of  current  research  is  directed  toward  the  implementation 
of  truly  multi-dimensional  upwind  schemes  in  which  the  upwind  biasing  is  determined  by  properties  of  the 
flow  rather  than  the  mesh  [69,  152,  109,  46,  171].  A  thorough  review  is  given  by  Pailliere  and  DeConinck  in 
reference  [144]. 

Residual  distribution  schemes  are  an  attractive  approach  for  triangulated  meshes.  In  these  the  residual 
defined  by  the  space  derivatives  is  evaluated  for  each  cell,  and  then  distributed  to  the  vertices  with  weights 
which  depend  on  the  direction  of  convection.  For  a  scalar  conservation  law  the  weights  can  be  chosen  to 
maintain  positivity  with  minimum  cross  diffusion  in  the  direction  normal  to  the  flow.  For  the  Euler  equations 
the  residual  can  be  linearized  by  assuming  that  the  parameter  vector  with  components  and  ^/pH 

varies  flnearly  over  the  cell.  Then 

dfjjw)  ^  .  dw 

dxj  ^  dxj 


where  the  Jacobian  matrices  Aj  =  ^  are  evaluated  with  Roe  averaging  of  the  values  of  to  at  the  vertices. 
Waves  in  the  direction  n  can  then  be  expressed  in  terms  of  the  eigenvectors  of  rijAj ,  and  a  positive  distribution 
scheme  is  used  for  waves  in  preferred  directions.  The  best  dioice  of  these  directions  is  the  subject  of  ongoing 


20 


research,  but  preliminary  results  indicate  the  possibility  of  achieving  high  resolution  of  shocks  and  contact 
discontinuities  which  are  not  aligned  with  mesh  lines  [144]. 

Hirsch  and  Van  Ransbeeck  adopt  an  alternative  approach  in  which  they  directly  construct  directional 
diffusive  terms  on  structured  meshes,  with  anti-diffusion  controlled  by  limiters  based  on  comparisons  of 
slopes  in  different  directions  [70].  They  also  show  promising  results  in  calculations  of  nozzles  with  multiply 
reflected  oblique  shocks. 

3.6.  Discretization  of  the  Viscous  Terms.  The  discretization  of  the  viscous  terms  of  the  Navier 
Stokes  equations  requires  an  approximation  to  the  velocity  derivatives  in  order  to  calculate  the  tensor 
aij,  equation  (3).  Then  the  viscous  terms  may  be  included  in  the  flux  balance  (4).  In  order  to  evaluate  the 
derivatives  one  may  apply  the  Gauss  formula  to  a  control  volume  V  with  the  boundary  S. 

f  =  [  UiUjds 

Jv  9xj  Js 

where  nj  is  the  outward  normal.  For  a  tetrahedral  or  hexahedral  cell  this  gives 


(38) 


dui  _  1 

dxj  vol 

faces 


Ui  rij  $ 


where  Ui  is  an  estimate  of  the  average  of  Ui  over  the  face.  If  u  varies  linearly  over  a  tetrahedral  cell  this  is 
exact.  Alternatively,  assuming  a  local  transformation  to  computational  coordinates  one  may  apply  the 
chain  rule 


(39) 


du 

dx 


dv' 

'2i' 

du 

dx 

M. 

dx 

Ico 

II 

M. 

Here  the  transformation  derivatives  ^  can  be  evaluated  by  the  same  finite  difference  formulas  as  the  velocity 
derivatives  In  this  case  ^  is  exact  if  u  is  a  linearly  varying  function. 

For  a  cell-centered  discretization  (figure  6a)  ^  is  needed  at  each  face.  The  simplest  procedure  is 
to  evaluate  ^  in  each  cell,  and  to  average  ^  between  the  two  cells  on  either  side  of  a  face  [95].  The 
resulting  discretization  does  not  have  a  compact  stencil,  and  supports  imdamped  oscillatory  modes.  In  a 
one  dimensional  calculation,  for  example,  0  would  be  discretized  as  ■  In  order  to  prodtice  a 

compact  stencil  may  be  estimated  from  a  control  volume  centered  on  each  face,  using  formulas  (38)  or 
(39)  [160].  This  is  computationally  expensive  because  the  number  of  faces  is  much  larger  than  the  number 
of  cells.  In  a  hexahedral  mesh  with  a  large  munber  of  vertices  the  munber  of  faces  approaches  three  times 
the  number  of  cells. 

This  motivates  the  introduction  of  dual  meshes  for  the  evaluation  of  the  velocity  derivatives  and  the  flux 
balance  as  sketched  in  figure  6.  The  figure  shows  both  cell-centered  and  cell- vertex  schemes.  The  dual  mesh 
connects  ceU  centers  of  the  primary  mesh.  If  there  is  a  kink  in  the  primary  mesh,  the  dual  cells  should  be 
formed  by  assembhng  contiguous  fractions  of  the  neighboring  primary  cells.  On  smooth  meshes  comparable 
results  are  obtained  by  either  of  these  formulations  [122,  123,  115].  If  the  mesh  has  a  kink  the  cell- vertex 
scheme  has  the  advantage  that  the  derivatives  are  calculated  in  the  interior  of  a  regular  cell,  with  no 
loss  of  accuracy. 

A  desirable  property  is  that  a  linearly  varying  velocity  distribution,  as  in  a  Couette  fiow ,  should  produce 
a  constant  stress  and  hence  an  exact  stress  balance.  This  property  is  not  necessarily  satisfied  in  general 
by  finite  difference  or  finite  volume  schemes  on  curvilinear  meshes.  The  characterization  fc-exact  has  been 
proposed  for  schemes  that  are  exact  for  polynomials  of  degree  k.  The  cell-vertex  finite  volume  scheme  is 
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centered  scheme.  aij 
evaluated  at  vertices  of 
the  primary  mesh 


Figure  6b:  Cell- 

vertex  scheme,  aij  evalu¬ 
ated  at  cell  centers  of  the 
primary  mesh 


Fig.  6.  Viscous  discretizations  for  cell- centered  and  cell-vertex  algorithms. 


linearly  exact  if  the  derivatives  are  evaluated  by  equation  (39),  since  then  is  exactly  evaluated  as  a 
constant,  leading  to  constant  viscous  stresses  aij,  and  an  exact  viscous  stress  balance.  This  remains  true 
when  there  is  a  kink  in  the  mesh,  because  the  summation  of  constant  stresses  over  the  faces  of  the  kinked 
control  volume  sketched  in  figure  6  still  yields  a  perfect  balance.  The  use  of  equation  (39)  to  evaluate  |^, 
however,  requires  the  additional  calculation  or  storage  of  the  nine  metric  quantities  in  each  cell,  whereas 
equation  (38)  can  be  evaluated  from  the  same  face  areas  that  are  used  for  the  flux  balance. 

In  the  case  of  an  unstructured  mesh,  the  weak  form  (6)  leads  to  a  natural  discretization  with  linear 
elements,  in  which  the  piecewise  linear  approximation  yields  a  constant  stress  in  each  ceU.  This  method 
yields  a  representation  which  is  globally  correct  when  averaged  over  the  cells,  as  is  proved  by  energy  estimates 
for  elliptic  problems  [19].  It  should  be  noted,  however,  that  it  yields  formulas  that  are  not  necessarily  locally 
consistent  with  the  differential  equations,  if  Taylor  series  expansions  are  substituted  for  the  solution  at  the 
vertices  appearing  in  the  local  stencil.  Figure  7  illustrates  the  discretization  of  the  Laplacian  u^x  + 
which  is  obtained  with  linear  elements.  It  shows  a  particular  triangulation  such  that  the  approximation  is 
locally  consistent  with  Uxx  +  Suyy.  Thus  the  use  of  an  irregular  triangulation  in  the  boundary  layer  may 
significantly  degrade  the  accuracy. 
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Fig.  7.  Example  of  discretization  Uxx  +  VLyy  on  a  triangxdar  mesh. 
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The  discretization  is  locally  equivalent  to  the  approxi- 
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Anisotropic  grids  are  needed  in  order  to  resolve  the  thin  boundary  layers  which  appear  in  viscous  flows 
at  high  Reynolds  numbers.  Otherwise  an  excessively  large  number  of  grid  cells  may  be  required.  The  use  of 
flat  tetrahedra  can  have  an  adverse  effect  on  both  the  accuracy  of  the  solution  and  the  rate  of  convergence 
to  a  steady  state.  This  has  motivated  the  use  of  hybrid  prismatic-tetrahedral  grids  in  which  prismatic  cells 
are  used  in  the  wall  regions  [145].  A  review  of  many  of  the  key  issues  in  the  design  of  flow  solvers  for 
imstructured  meshes  is  given  by  Venkatakrishnan  [194]. 

3.7.  Time  Stepping  Schemes.  If  the  space  discretization  procedure  is  implemented  separately,  it 
leads  to  a  set  of  coupled  ordinary  differential  equations,  which  can  be  written  in  the  form 

(40) 

where  w  is  the  vector  of  the  flow  variables  at  the  mesh  points,  and  R.(w)  is  the  vector  of  the  residuals, 
consisting  of  the  flux  balances  defined  by  the  space  discretization  scheme,  together  with  the  added  dissipative 
terms.  If  the  objective  is  simply  to  reach  the  steady  state  and  details  of  the  transient  solution  are  immaterial, 
the  time-stepping  scheme  may  be  designed  solely  to  maximize  the  rate  of  convergence.  The  first  decision  that 
must  be  made  is  whether  to  use  an  explicit  scheme,  in  which  the  space  derivatives  are  calculated  fi-om  known 
values  of  the  flow  variables  at  the  beginning  of  the  time  step,  or  an  implicit  scheme,  in  which  the  formulas  for 
the  space  derivatives  include  as  yet  unknown  values  of  the  flow  variables  at  the  end  of  the  time  step,  leading 
to  the  need  to  solve  coupled  equations  for  the  new  values.  The  permissible  time  step  for  an  explicit  scheme 
is  limited  by  the  Courant-Priedrichs-Lewy  (CFL)  condition,  which  states  that  a  difference  scheme  cannot  be 
a  convergent  and  stable  approximation  unless  its  domain  of  dependence  contains  the  domain  of  dependence 
of  the  corresponding  differential  equation.  One  can  anticipate  that  implicit  schemes  will  yield  convergence 
in  a  smaller  number  of  time  steps,  because  the  time  step  is  no  longer  constrained  by  the  CFL  condition. 
Implicit  schemes  will  be  eflficient,  however,  only  if  the  decrease  in  the  number  of  time  steps  outweighs  the 
increase  in  the  computational  effort  per  time  step  consequent  upon  the  need  to  solve  coupled  equations.  The 
prototype  implicit  scheme  can  be  formulated  by  estimating  at  f  +  //Af  as  a  linear  combination  of  R.(w”) 
and  R(w"+^).  The  resulting  equation 

w"+i  =  w”  -  Af  {(1  -  /i)  R  (w")  -I-  fiR  (w"'+^) } 

can  be  linearized  as 

+  nAt^^Sw  +  AtR  (w^)  =  0. 

If  one  sets  =  1  and  lets  Af  — >■  oo  this  reduces  to  the  Newton  iteration  ,  which  has  been  successfully  used 
in  two-dimensional  calculations  [192,  59].  In  the  three-dimensional  case  with,  say,  an  N  x  N  x  N  mesh, 
the  bandwidth  of  the  matrix  that  must  be  inverted  is  of  order  N^.  Direct  inversion  requires  a  number  of 
operations  proportional  to  the  number  of  unknowns  multiplied  by  the  square  of  the  bandwidth  of  the  order 
of  N"^.  This  is  prohibitive,  and  forces  recourse  to  either  an  approximate  factorization  method  or  an  iterative 
solution  method. 

Alternating  direction  methods,  which  introduce  factors  corresponding  to  each  coordinate,  are  widely  used 
for  structured  meshes  [21,  154].  They  cannot  be  implemented  on  unstructured  tetrahedral  meshes  that  do 
not  contain  identifiable  mesh  directions,  although  other  decompositions  are  possible  [66,  116].  If  one  chooses 
to  adopt  the  iterative  solution  technique,  the  principal  alternatives  are  variants  of  the  Gauss-Seidel  and 
Jacobi  methods.  A  symmetric  Gauss-Seidel  method  with  one  iteration  per  time  step  is  essentially  equivalent 
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to  an  approximate  lower-iipper  (LU)  factorization  of  the  implicit  scheme  [94,  137,  36,  205].  On  the  other 
hand,  the  Jacobi  method  with  a  fixed  number  of  iterations  per  time  step  reduces  to  a  multistage  explicit 
scheme,  belonging  to  the  general  class  of  Runge-Kutta  schemes  [38].  Schemes  of  this  type  have  proved  very 
effective  for  wide  variety  of  problems,  and  they  have  the  advantage  that  they  can  be  applied  equally  easily 
on  both  structured  and  unstnictured  meshes  [93,  78,  80,  161]. 

If  one  reduces  the  linear  model  problem  corresponding  to  (40)  to  an  ordinary  differential  equation  by 
substituting  a  Fourier  mode  w  =  the  resTilting  Fourier  symbol  has  an  imaginary  part  proportional  to 
the  wave  speed,  and  a  negative  real  part  proportional  to  the  diffusion.  Thus  the  time  stepping  scheme  should 
have  a  stability  region  which  contains  a  substantial  interval  of  the  negative  real  axis,  as  well  as  an  interval 
along  the  imaginary  axis.  To  achieve  this  it  pays  to  treat  the  convective  and  dissipative  terms  in  a  distinct 
fashion.  Thus  the  residual  is  spUt  as 

R{w)  =  Q{w)  +  D{w), 

where  Q{w)  is  the  convective  part  and  D{w)  the  dissipative  part.  Denote  the  time  level  nAt  by  a  superscript 
n.  Then  the  multistage  time  stepping  scheme  is  formulated  as 

^(n+1,0)  ^ 

^(n+l.fc)  =  _  akAt 


where  the  superscript  k  denotes  the  fc-th  stage,  Om  =  R  and 

=  Q  (m”) ,  =  D  (w”) 


+  (1  - 

The  coefficients  ak  are  chosen  to  maximize  the  stability  interval  along  the  imaginary  axis,  and  the  coefficients 
Pk  are  chosen  to  increase  the  stability  interval  along  the  negative  real  axis. 

These  schemes  do  not  fall  within  the  standard  framework  of  Runge-Kutta  schemes,  and  they  have  much 
larger  stability  regions  [80] .  Two  schemes  which  have  been  found  to  be  particularly  effective  are  tabulated 
below.  The  first  is  a  four-stage  scheme  with  two  evaluations  of  dissipation.  Its  coefficients  are 
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The  second  is  a  five-stage  scheme  with  three  evaluations  of  dissipation.  Its  coefficients  are 
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3.8.  Multigrid  Methods. 

3.8.1.  Acceleration  of  Steady  Flow  Calculations.  Radical  improvements  in  the  rate  of  convergence 
to  a  steady  state  can  be  realized  by  the  miiltigrid  time-stepping  technique.  The  concept  of  acceleration  by  the 
introduction  of  multiple  grids  was  first  proposed  by  Fedorenko  [57],  There  is  by  now  a  fairly  well-developed 
theory  of  multigrid  methods  for  elliptic  equations  based  on  the  concept  that  the  updating  scheme  acts  as  a 
smoothing  operator  on  each  grid  [28,  62].  This  theory  does  not  hold  for  hyperbolic  systems.  Nevertheless, 
it  seems  that  it  ought  to  be  possible  to  accelerate  the  evolution  of  a  hyperbolic  system  to  a  steady  state 
by  using  large  time  steps  on  coarse  grids  so  that  disturbances  will  be  more  rapidly  expelled  through  the 
outer  boundary.  Various  multigrid  time-stepping  schemes  designed  to  take  advantage  of  this  effect  have 
been  proposed  [136,  77,  64,  81,  34,  9,  67] 

One  can  devise  a  multigrid  scheme  using  a  sequence  of  independently  generated  coarser  meshes  by 
ehminating  alternate  points  in  each  coordinate  direction.  In  order  to  give  a  precise  description  of  the 
multigrid  scheme,  subscripts  may  be  used  to  indicate  the  grid.  Several  transfer  operations  need  to  be 
defined.  First  the  solution  vector  on  grid  k  must  be  initiahzed  as 


where  Wk-i  is  the  current  value  on  grid  k  —  I,  and  Tk,k-\  is  a  transfer  operator.  Next  it  is  necessary  to 
transfer  a  residual  forcing  function  such  that  the  solution  grid  k  is  driven  by  the  residuals  calculated  on  grid 
A;  —  1.  This  can  be  accomplished  by  setting 

Pk  —  Qk,k-lRk-l  (wfe-l)  —  Rk  , 

where  Qk,k-i  is  another  transfer  operator.  Then  Rk{wk)  is  replaced  by  Rk{wk)  +  P/t  in  the  time-  stepping 
scheme.  Thus,  the  multistage  scheme  is  reformulated  as 

-  aiAtk  +  Pfej 


-  Olq+lAtk 


Rk^+Pk  ■ 


The  result  then  provides  the  initial  data  for  grid  k  +  l.  Finally,  the  accmnulated  correction  on  grid 

k  has  to  be  transferred  back  to  grid  k-l  with  the  aid  of  an  interpolation  operator  Ik-i,k-  With  properly 
optimized  coefficients  multistage  time-stepping  schemes  can  be  very  efficient  drivers  of  the  multigrid  process. 
A  W-cycle  of  the  type  illustrated  in  Figure  8  proves  to  be  a  particularly  effective  strategy  for  managing  the 
work  split  between  the  meshes.  In  a  three-dimensional  case  the  number  of  cells  is  reduced  by  a  factor  of 
eight  on  each  coarser  grid.  On  examination  of  the  figure,  it  can  therefore  be  seen  that  the  work  measured  in 
imits  corresponding  to  a  step  on  the  fine  grid  is  of  the  order  of 


1  -f  2/8  -I-  4/64  +  ...<  4/3, 


and  consequently  the  very  large  effective  time  step  of  the  complete  cycle  costs  only  slightly  more  than  a 
single  time  step  in  the  fine  grid. 

This  procedmre  has  proved  extremely  successful  for  the  solution  of  the  inviscid  Euler  equations,  but  less 
effective  in  calculations  of  turbulent  viscous  flows  at  high  Reynolds  numbers  using  the  Reynolds  averaged 
Navier  Stokes  equations.  These  require  highly  anisotropic  grids  with  very  fine  mesh  intervals  normal  to  the 
waU  to  resolve  the  boundary  layers.  While  simple  multigrid  methods  stiU  yield  fast  initial  convergence,  they 
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Fig.  8.  Multigrid  W -cycle  for  managing  the  grid  calculation.  E,  evaluate  the  change  in  the  flow  for  one  step;  T,  transfer 
the  data  without  updating  the  solution. 


tend  to  slow  down  as  the  calculation  proceeds  to  a  low  asymptotic  rate.  This  has  motivated  the  introduction 
of  semi-coarsening  and  directional  coarsening  methods  [132,  133,  3,  4,  5,  148,  149]. 

The  multigrid  method  can  be  apphed  on  unstnictmed  meshes  by  interpolating  between  a  sequence 
of  separately  generated  meshes  with  progressively  increasing  cell  sizes  [91,  126,  127,  146).  It  is  not  easy 
to  generate  very  coarse  meshes  for  complex  configurations.  An  alternative  approach,  which  removes  this 
difficulty,  is  to  automatically  generate  successively  coarser  meshes  by  agglomerating  control  volumes  or  by 
collapsing  edges.  This  approach  yields  comparable  rates  of  convergence  and  has  proved  to  be  quite  robust 
[101,  102,  128,  42]. 

3.8.2.  Multigrid  Implicit  Schemes  for  Unsteady  Flow.  Time  dependent  calculations  are  needed 
for  a  number  of  important  applications,  such  as  flutter  analysis,  or  the  analysis  of  the  flow  past  a  helicopter 
rotor,  in  which  the  stabihty  limit  of  an  explicit  scheme  forces  the  use  of  much  smaller  time  steps  than  would 
be  needed  for  an  accurate  simulation.  In  this  situation  a  multigrid  explicit  scheme  can  be  used  in  an  inner 
iteration  to  solve  the  equations  of  a  fully  impUcit  time  stepping  scheme  [84]. 

Suppose  that  (40)  is  approximated  as 

+  R{w^+'^)  =  0. 


Here  Dt  is  a  order  accurate  backward  difference  operator  of  the  form 

q=l  ^ 


where 
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Applied  to  the  linear  differential  equation 

dw 

“  aw 
dt 

the  schemes  with  k  =  1,2  are  stable  for  all  o^At  in  the  left  half  plane  (A-stable),  Dahlquist  has  shown  that 
A-stable  linear  multi-step  schemes  are  at  best  second  order  accurate  [44].  Gear  however,  has  shown  that  the 
schemes  with  fc  <  6  are  stiffly  stable  [58],  and  one  of  the  higher  order  schemes  may  offer  a  better  compromise 
between  accuracy  and  stability,  depending  on  the  application. 

Equation  (40)  is  now  treated  as  a  modified  steady  state  problem  to  be  solved  by  a  multigrid  scheme 
using  variable  local  time  steps  in  a  fictitious  time  t* ,  For  example,  in  the  case  k  —  2  one  solves 


where 


2At 


W 


1 


and  the  last  two  terms  are  treated  as  fixed  source  terms.  The  first  term  shifts  the  Fourier  symbol  of  the 
equivalent  model  problem  to  the  left  in  the  complex  plane.  While  this  promotes  stability,  it  may  also  require 
a  limit  to  be  imposed  on  the  magnitude  of  the  local  time  step  At*  relative  to  that  of  the  implicit  time  step 
At.  This  may  be  relieved  by  a  point-implicit  modification  of  the  multi-stage  scheme  [129].  In  the  case  of 
problems  with  moving  boundaries  the  equations  must  be  modified  to  allow  for  movement  and  deformation 
of  the  mesh. 

This  method  has  proved  effective  for  the  calculation  of  unsteady  flows  that  might  be  associated  with 
wing  fiutter  [6,  7]  and  also  in  the  calculation  of  unsteady  incompressible  flows  [22].  It  has  the  advantage  that 
it  can  be  added  as  an  option  to  a  computer  program  which  uses  an  explicit  multigrid  scheme,  allowing  it  to  be 
used  for  the  efficient  calculation  of  both  steady  and  imsteady  flows.  A  similar  approach  has  been  successfully 
adopted  for  unsteady  flow  simulations  on  unstructured  grids  by  Venkatakrishnan  and  Mavriplis  [195]. 


3.9.  Preconditioning.  Another  way  to  improve  the  rate  of  convergence  to  a  steady  state  is  to  mul¬ 
tiply  the  space  derivatives  in  equation  (1)  by  a  preconditioning  matrix  P  which  is  designed  to  equalize  the 
eigenvalues,  so  that  aU  the  waves  can  be  advanced  with  optimal  time  steps.  A  symmetric  preconditioner 
which  equalizes  the  eigenvalues  has  been  proposed  by  Van  Leer  [110].  When  the  equations  are  written  in 
stream-aligned  coordinates  this  has  the  form 


where 


P  = 


0  0  0 

-fM  ^  +  1  0  0  0 

0  0  r  0  0 

0  0  0  T  0 

0  0  0  0  1 


^  =  T  =  \/l  —  if  Af  <  1 

if  M>1 


Turkel  has  proposed  an  as5Tnmetric  preconditioner  which  has  also  proved  effective,  particularly  for  flow  at 
low  Mach  numbers  [191].  The  use  of  these  preconditioners  can  lead  to  instability  at  stagnation  points  where 
there  is  a  zero  eigenvalue  which  cannot  be  equalized  with  the  eigenvalues  ±c. 
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The  preconditioners  of  Van  Leer  and  Turkel  do  not  take  account  of  the  effect  of  differences  in  the 
mesh  intervals  in  the  different  coordinate  directions.  The  need  to  resolve  the  boundary  layer  generally 
compels  the  introduction  of  mesh  cells  with  very  high  aspect  ratios  near  the  boundary,  and  these  can  lead 
to  a  severe  reduction  in  the  rate  of  convergence  to  a  steady  state.  AUmaras  has  analyzed  explicit  and 
implicit  Jacobi-based  preconditioners  which  include  the  influence  of  the  mesh  intervals  [3,  4,  5],  Using 
a  block-Jacobi  preconditioner  with  coarsening  only  in  the  direction  normal  to  the  wall,  Pierce  has  recently 
obtained  impressive  results  on  viscous  flows  with  high  aspect  ratio  grids  [148, 149].  Mavriplis  has  successfully 
combined  block  preconditioners  with  line  solvers  to  accelerate  the  convergence  of  viscous  flow  solutions  on 
highly  stretched  imstructured  grids  [125]. 

An  alternative  approach  has  recently  been  proposed  by  Ta’asan  [186],  in  which  the  equations  are  written 
in  a  canonical  form  which  separates  the  equations  describing  acoustic  waves  from  those  describing  convection. 
In  terms  of  the  velocity  components  u,v  and  the  vorticity  w,  temperature  T,  entropy  s  and  total  enthalpy 
H,  the  equations  describing  steady  two-dimensional  flow  can  be  written  as 

■  T»i 

__d_ 

dy 

0 
0 
0 

where 


and 


Here  the  first  three  equations  describe  an  elliptic  system  if  the  flow  is  subsonic,  while  the  remaining  equations 
are  convective.  Now  separately  optimized  multigrid  procedures  are  used  to  solve  the  two  sets  of  equations, 
which  are  essentially  decoupled.  An  alternative  approach  to  the  optimal  splitting  of  the  flow  equations  into 
convective  and  acoustic  parts  has  been  developed  by  Sidilkover  [172,  162]. 

3,10.  High  Order  Schemes  and  Mesh  Refinement.  The  need  both  to  improve  the  accuracy  of 
computational  simulations,  and  to  assure  known  levels  of  accuracy  is  the  focus  of  ongoing  research.  The 
main  routes  to  improving  the  accuracy  are  to  increase  the  order  of  the  discrete  scheme,  and  reduce  the  mesh 
interval.  High  order  difference  methods  are  most  easily  implemented  on  Cartesian,  or  at  least  extremely 
smooth  grids.  The  expansion  of  the  stencil  as  the  order  is  increased  leads  to  the  need  for  complex  boundary 
conditions.  Compact  schemes  keep  the  stencil  as  small  as  possible  [157,  112,  33].  On  simple  domains, 
spectral  methods  are  particularly  effective,  especially  in  the  case  of  periodic  boundary  conditions,  and  can 
be  used  to  produce  exponentially  fast  convergence  of  the  error  as  the  mesh  interval  is  decreased  [139,  32].  A 
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compromise  is  to  divide  the  field  into  subdomains  and  introduce  high  order  elements.  This  approach  is  used 
in  the  spectral  element  method  [100]. 

High  order  difference  schemes  and  spectral  methods  have  proven  particularly  useful  in  direct  Navier- 
Stokes  simulations  of  transient  and  turbulent  flows.  High  order  niethods  are  also  beneficial  in  computational 
aero-acoustics,  where  it  is  desired  to  track  waves  over  long  distances  with  minimum  error.  If  the  flow  contains 
shock  waves  or  contact  discontinuities,  the  ENO  method  may  be  used  to  construct  high  order  non-oscillatory 
schemes. 

In  multi-dimensional  flow  simulations,  global  reduction  of  the  mesh  interval  can  be  prohibitively  expen¬ 
sive,  motivating  the  use  of  adaptive  mesh  refinement  procedures  which  reduce  the  local  mesh  width  h  if  there 
is  an  indication  that  the  error  is  too  large  [25,  45,  117,  71,  155,  111,  147].  In  such  h-reflnement  methods, 
simple  error  indicators  such  as  local  solution  gradients  may  be  used.  Alternatively,  the  discretization  error 
may  be  estimated  by  comparing  quantities  calculated  with  two  mesh  widths,  say  on  the  current  mesh  and  a 
coarser  mesh  with  double  the  mesh  interval.  Procedures  of  this  kind  may  also  be  used  to  provide  a  posteriori 
estimates  of  the  error  once  the  calculation  is  completed. 

This  kind  of  local  adaptive  control  can  also  be  applied  to  the  local  order  of  a  finite  element  method 
to  produce  a  p-refinement  method,  where  p  represents  the  order  of  the  polynomial  basis  functions.  Finally, 
both  /i-  and  p-  refinement  can  be  combined  to  produce  an  h-p  method  in  which  h  and  p  are  locally  optimized 
to  yield  a  solution  with  minimum  error  [138].  Such  methods  can  achieve  exponentially  fast  convergence,  and 
are  weU  established  in  computational  solid  mechanics. 

4.  Current  Status  of  Numerical  Simulation*  This  section  presents  some  representative  numerical 
results  which  confirm  the  properties  of  the  algorithms  which  have  been  reviewed  in  the  last  section.  These 
have  been  drawn  from  the  work  of  the  author  and  his  associates.  They  also  illustrate  the  kind  of  calculation 
which  can  be  performed  in  an  industrial  environment,  where  rapid  turn  around  is  important  to  allow  the 
quick  assessment  of  design  changes,  and  computational  costs  must  be  limited. 

4*1.  One  dimensional  shock*  In  order  to  verify  the  discrete  structure  of  stationary  shocks,  calcu¬ 
lations  were  performed  for  a  one  dimensional  problem  with  initial  data  containing  left  and  right  states 
compatible  with  the  Rankine  Hugoniot  conditions.  An  intermediate  state  consisting  of  the  arithmetic  av¬ 
erage  of  the  left  and  right  states  was  introduced  at  a  single  ceU  in  the  center  of  the  domain.  With  this 
intermediate  state  the  system  is  not  in  equilibrium,  and  the  time  dependent  equations  were  solved  to  find 
an  equilibrium  solution  with  a  stationary  shock  wave  separating  the  left  and  right  states.  Table  1  shows  the 
result  for  a  shock  wave  at  Mach  20.  This  calculation  used  the  H-CUSP  scheme,  which  allows  a  solution  with 
constant  stagnation  enthalpy.  The  SLIP-JST  construction  was  used  with  the  limiter  defined  by  equation 
(23),  and  g  =  3.  The  table  shows  the  values  of  p,  n,  if,  p,  M  and  the  entropy  5  =  log^  —  log 
A  perfect  one  point  shock  structure  is  displayed.  The  entropy  is  zero  to  4  decimal  places  upstream  of  the 
shock,  exhibits  a  slight  excursion  at  the  interior  point,  and  is  constant  to  4  decimal  places  downstream  of  the 
shock.  It  may  be  noted  that  the  mass,  momentum  and  energy  of  the  initial  data  are  not  compatible  with  the 
final  equilibrium  state.  According  to  conservation  arguements  the  total  mass,  momentum  and  energy  must 
remain  constant  if  the  outflow  flux  /h  remains  equal  to  the  inflow  flux  /l-  Therefore  fji  must  be  allowed 
to  vary  according  to  an  appropriate  outflow  boundary  condition  to  allow  the  total  mass,  momentum  and 
energy  to  be  adjusted  to  values  compatible  with  equilibrium. 

4.2*  Euler  Calculations  for  Airfoils  and  Wings.  The  results  of  transonic  flow  calculations  for  two 
well  known  airfoils,  the  RAE  2822  and  the  NACA  0012,  are  presented  in  figures  (21-23).  The  H-CUSP 
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P 

H 

P 

M 

19 

1.0000 

283.5000 

1.0000 

20.0000 

20 

1.0000 

283.5000 

1.0000 

20.0000 

21 

1.0000 

283.5000 

1.0000 

20.0000 

22 

4.1924 

283.4960 

307.4467 

0.7229 

23 

5.9259 

283.4960 

466.4889 

0.3804 

24 

5.9259 

283.4960 

466.4889 

0.3804 

25 

5.9259 

283.4960 

466.4889 

0.3804 

Table  1 

Shock  Wave  at  Mach  20 


scheme  was  again  used  with  the  SLIP-JST  construction.  The  limiter  defined  by  equation  (23)  was  used  with 
q  =  3.  The  5  stage  time  stepping  scheme  (42)  was  augmented  by  the  multigrid  scheme  described  in  section 
4.2  to  accelerate  convergence  to  a  steady  state.  The  equations  were  discretized  on  meshes  with  0-topology 
extending  out  to  a  radius  of  about  100  chords.  In  each  case  the  calculations  were  performed  on  a  sequence 
of  successively  finer  meshes  from  40x8  to  320x64  cells,  while  the  multigrid  cycles  on  each  of  these  meshes 
descended  to  a  coarsest  mesh  of  10x2  cells.  Figures  21-23  show  the  final  results  on  320x64  meshes  for  the 
RAE  2822  airfoil  at  Mach  .75  and  3°  angle  of  attack,  and  for  the  NACA  0012  airfoil  at  Mach  .8  and  1.25° 
angle  of  attack,  and  also  at  Mach  .85  and  1°  angle  of  attack.  In  the  pressure  distributions  the  pressure 
coefficient  is  plotted  with  the  negative  (suction)  pressures  upward,  so  that  the  upper  curve 

represents  the  flow  over  the  upper  side  of  a  lifting  airfoil.  The  convergence  histories  show  the  mean  rate 
of  change  of  the  density,  and  also  the  total  number  of  supersonic  points  in  the  flow  field,  which  provides 
a  useful  measure  of  the  global  convergence  of  transonic  flow  calculations  such  as  these.  In  each  case  the 
convergence  history  is  shown  for  100  cycles,  while  the  pressure  distribution  is  displayed  after  a  sufficient 
number  of  cycles  for  its  convergence.  The  pressure  distribution  of  the  RAE  2822  airfoil  converged  in  only  25 
cycles.  Convergence  was  slower  for  the  NACA  0012  airfoil.  In  the  case  of  flow  at  Mach  .8  and  1.25°  angle  of 
attack,  additional  cycles  were  needed  to  damp  out  a  wave  downstream  of  the  weak  shock  wave  on  the  lower 
surface. 

As  a  further  check  on  accuracy  the  drag  coefficient  should  be  zero  in  subsonic  flow,  or  in  shock  free 
transonic  flow.  Table  2  shows  the  computed  drag  coefficient  on  a  sequence  of  three  meshes  for  three  examples. 
The  first  two  are  subsonic  flows  over  the  RAE  2822  and  NACA  0012  airfoils  at  Mach  .5  and  3°  angle  of 
attack.  The  third  is  the  flow  over  the  shock  free  Korn  airfoil  at  its  design  point  of  Mach  .75  and  0°  angle  of 
attack.  In  all  three  cases  the  drag  coefficient  is  calculated  to  be  zero  to  four  digits  on  a  160x32  mesh. 


Mesh 

RAE  2822 

Mach  .50 

a  3° 

NACA  0012 

Mach  .50 

a  3° 

Korn  Airfoil 

Mach  .75 

40x8 

.0062 

.0047 

.0098 

80x16 

.0013 

.0008 

.0017 

160x32 

.0000 

.0000 

.0000 

Table  2 

Drag  Coefficient  on  a  sequence  of  meshes 
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As  a  further  test  of  the  performance  of  the  H-CUSP  scheme,  the  flow  past  the  ONERA  M6  wing  was 
calculated  on  a  mesh  with  C-H  topology  and  192x32x48  =  294912  cells.  Figure  24  shows  the  result  at  Mach 
.84  and  3.06°  angle  of  attack.  This  again  verifies  the  non-oscillatory  character  of  the  solution,  and  the  sharp 
resolution  of  shock  waves.  In  this  case  50  cycles  were  sufficient  for  convergence  of  the  pressure  distributions. 

Figure  9  shows  a  calculation  of  the  Northrop  YF23  by  R.J.  Busch,  Jr.,  who  used  the  author’s  FL057 
code  to  solve  the  Euler  equations  [30].  Although  an  inviscid  model  of  the  flow  was  used,  it  can  be  seen  that 
the  simulations  are  in  good  agreement  with  wind  tunnel  measurements  both  at  Mach  .90,  with  angles  of 
attack  of  0,  8  and  16  degrees,  and  at  Mach  1.5  with  angles  of  attack  of  0,  4  and  8  degrees.  At  a  high  angle  of 
attack  the  flow  separates  from  the  leading  edge,  and  this  example  shows  that  in  situations  where  the  point 
of  separation  is  fixed,  an  inviscid  model  may  still  produce  a  useful  prediction.  Thus  valuable  information  for 
the  aerodynamic  design  could  be  obtained  with  a  relatively  inexpensive  computational  model. 


Fig.  9.  Comparison  of  Experimented  and  Computed  Drag  Rise  Curve  for  the  YF-ZS  (Supplied  by  R.  J.  Bush  Jr.) 

The  next  figures  show  the  results  of  calculations  using  the  AIRPLANE  code  developed  by  T.  J.  Baker 
and  the  author,  to  solve  the  Euler  equations  on  an  unstructured  mesh.  This  provides  the  flexibility  to 
treat  arbitrarily  complex  configurations  without  the  need  to  spend  months  developing  an  acceptable  mesh. 
Figures  10  and  11  show  calculations  for  supersonic  transport  configurations  which  were  performed  by  Susan 
Cliff.  The  agreement  with  experimental  data  is  quite  good,  and  it  has  also  been  possible  to  predict  the 
sonic  boom  signature  [39].  Figure  12  shows  an  Euler  calculation  for  the  McDonnell  Douglas  MDll  with  flow 
through  the  engine  nacelles,  using  348407  mesh  points  of  2100466  tetrahedra.  This  calculation  takes  4  hours 
on  an  IBM  590  workstation.  A  parallel  version  of  the  code  has  been  developed  in  collaboration  with  W.S. 
Chen,  and  the  same  calculation  can  be  performed  in  16  minutes  using  16  processors  of  an  IBM  SP2.  The 
parallel  speed-up  for  the  MDll  is  shown  in  table  3. 
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Fig.  10.  Gomparison  of  Experimental  and  Calculated  Results  for  a  HSCT  Configuration 


'  '  ' 


Fig.  11.  Pressure  Contours  and  Sonic  Boom  on  a  Representative  HSCT  Configuration 


4.3.  Viscous  Flow  Cedculations.  The  next  figures  show  viscous  simulations  based  on  the  solution  of 
the  Reynolds  averaged  Navier  Stokes  equations  with  turbulence  models.  Figure  13  shows  a  two-dimensional 
calculation  for  the  RAE  2822  airfoil  by  L.  Martinelli.  The  vertical  axis  represents  the  negative  pressure 
coefficient,  and  there  is  a  shock  wave  half  way  along  the  upper  surface.  This  example  confirms  that  in  the 
absence  of  significant  shock  induced  separation,  simulations  performed  on  a  sufficiently  fine  mesh  (512  x  64) 
can  produce  excellent  agreement  with  experimental  data.  Figure  20  shows  a  simulation  of  the  McDonneU- 
Douglas  F18  performed  by  R.M.  Cummings,  Y.M.  Rizk,  L.B.  Schiff  and  N.M.  Chaderjian  at  NASA  Ames 
[43].  They  used  a  multiblock  mesh  with  about  900000  mesh  points.  While  this  is  probably  not  enough  for 
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No.  of  Nodes 

Seconds/ Cycle 

Speedup 

1 

36.03 

1.00 

2 

18.11 

1.99 

4 

9.11 

3.96 

8 

4.66 

7.73 

16 

2.39 

15.08 

Table  3 

AIRPLANE  Parallel  Performance  on  the  SP2,  MD-11  Model 


an  accurate  quantitative  prediction,  the  agreement  with  both  the  experimental  data  and  the  visualization 
are  quite  good. 


Fig.  12.  Computed  Pressure  Field  for  a  McDonnell  Douglas  MDll 

Figure  14  shows  an  unsteady  flow  calculation  for  a  pitching  airfoil  performed  by  J.  Alonso  using  the 
code  UFL082,  which  he  jointly  developed  with  L.  Martinelli  and  the  author  [7].  This  uses  the  multigrid 
implicit  scheme  described  in  Section  3.7.2  which  allows  the  number  of  time  steps  to  be  reduced  from  several 
thousand  to  36  per  pitching  cycle.  The  agreement  with  experimental  data  is  quite  good. 

4.4.  Ship  Wave  Resistance  calculations.  Figures  51-53  show  the  results  of  an  application  of  the 
same  multigrid  finite  volume  techniques  to  the  calculation  of  the  flow  past  a  naval  frigate,  using  a  code  which 
was  developed  by  J.  Farmer,  L.  Martinelli  and  the  author  [56].  The  mesh  was  adjusted  during  the  course 
of  the  calculation  to  conform  to  the  free  surface  in  order  to  satisfy  the  exact  non-linear  boundary  condition, 
while  artificial  compressibility  was  used  to  reat  the  incompressible  flow  equations. 

5.  Aerodynamic  Shape  Optimization. 

5.1.  Optimization  and  Design.  Traditionally  the  process  of  selecting  design  variations  has  been 
carried  out  by  trial  and  error,  relying  on  the  intuition  and  experience  of  the  designer.  With  currently 
available  equipment  the  turn  around  for  numerical  simulations  is  becoming  so  rapid  that  it  is  feasible  to 
examine  an  extremely  large  number  of  variations.  It  is  not  at  all  hkely  that  repeated  trials  in  an  interactive 
design  and  analysis  procedure  can  lead  to  a  truly  optimum  design.  In  order  to  take  full  advantage  of  the 
possibility  of  examining  a  large  design  space  the  numerical  simulations  need  to  be  combined  with  automatic 
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Fig.  13.  Two-Dimensional  Turbulent  Viscous  Calculation  (by  Luigi  Martinelli) 


search  and  optimization  procedures.  This  can  lead  to  automatic  design  methods  which  will  fully  realize  the 
potential  improvements  in  aerod3mamic  efficiency. 

The  simplest  approach  to  optimization  is  to  define  the  geometry  through  a  set  of  design  parameters, 
which  may,  for  example,  be  the  weights  Oi  applied  to  a  set  of  shape  functions  bi{x)  so  that, the  shape  is 
represented  as 


/(^)  =^o!ibi(x). 

Then  a  cost  function  I  is  selected  which  might,  for  example,  be  the  drag  coefficient  at  a  final  lift  coefficient, 
and  I  is  regarded  as  a  function  of  the  parameters  .  The  sensitivities  ^  may  now  be  estimated  by  making 
a  small  variation  Sai  in  each  design  parameter  in  turn  and  recalculating  the  flow  to  obtain  the  change  in  J. 
Then 

dl  ^  I{ai  +  Scxj)  -  /(oj) 
dcxi  Sai 

The  gradient  vector  may  now  be  used  to  determine  a  direction  of  improvement.  The  simplest  procedure 
is  to  make  a  step  in  the  negative  gradient  direction  by  setting 

^n+l 


so  that  to  first  order 


1  +  51  =  1- 


dF  ^  ^  .dFdi 

da^°‘  ^  ^  da  da' 
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FiCt.  14.  Mach  Number  Contours.  Pitching  Airfoil  Case.  Re  =  1.0  x  10^,  Moo  —  0.796,  Kc  —  0.202. 


More  sophisticated  search  procedures  may  be  used  such  as  quasi-Newton  methods,  which  attempt  to  estimate 
the  second  derivative  of  the  cost  function  from  changes  in  the  gradient  in  successive  optimization 

steps.  These  methods  also  generally  introduce  line  searches  to  find  the  minimum  in  the  search  direction 
which  is  defined  at  each  step.  The  main  disadvantage  of  this  approach  is  the  need  for  a  number  of  flow 
calculations  proportional  to  the  number  of  design  variables  to  estimate  the  gradient.  The  computational 
costs  can  thus  become  prohibitive  as  the  number  of  design  variables  is  increased. 

An  alternative  approach  is  to  cast  the  design  problem  as  a  search  for  the  shape  that  will  generate  the 
desired  pressure  distribution.  This  approach  recognizes  that  the  designer  usually  has  an  idea  of  the  the 
kind  of  pressure  distribution  that  will  lead  to  the  desired  performance.  Thus,  it  is  useful  to  consider  the 
inverse  problem  of  calculating  the  shape  that  will  lead  to  a  given  pressure  distribution.  The  method  has  the 
advantage  that  only  one  flow  solution  is  required  to  obtain  the  desired  design.  Unfortunately,  a  physically 
realizable  shape  may  not  necessarily  exist,  unless  the  pressure  distribution  satisfies  certain  constraints.  Thus 
the  problem  must  be  very  carefully  formulated;  otherwise  it  may  be  ill  posed. 

The  difficulty  that  the  target  pressure  may  be  unattainable  may  be  circumvented  by  treating  the  inverse 
problem  as  a  special  case  of  the  optimization  problem,  with  a  cost  function  which  measures  the  error  in  the 
solution  of  the  inverse  problem.  For  example,  if  pd  is  the  desired  surface  pressure,  one  may  take  the  cost 
ftinction  to  be  an  integral  over  the  the  body  surface  of  the  square  of  the  pressure  error, 

I  j^ip-Pd)'^dB, 

or  possibly  a  more  general  Sobolev  norm  of  the  pressure  error.  This  has  the  advantage  of  converting  a 
possibly  ill  posed  problem  into  a  well  posed  one.  It  has  the  disadvantage  that  it  incurs  the  computational 
costs  associated  with  optimization  procedures. 

5.2.  Application  of  Control  Theory.  In  order  to  rediice  the  computational  costs,  it  turns  out  that 
there  are  advantages  in  formulating  both  the  inverse  problem  and  more  general  aerodynamic  problems 


Fig.  15,  Contours  of  Surf  ace  Wave  Elevation  for  a  Combatant  Ship 


Fig.  16.  Contours  of  Surface  Wave  Elevation  Near  the  Transom  Stem 


within  the  framework  of  the  mathemaitical  theory  for  the  control  of  systems  governed  by  partial  differential 
equations  [113],  A  wing,  for  example,  is  a  device  to  produce  lift  by  controlling  the  flow,  and  its  design 
can  be  regarded  as  a  problem  in  the  optimal  control  of  the  flow  equations  by  variation  of  the  shape  of  the 
boundary.  If  the  boundary  shape  is  regarded  as  arbitrary  within  some  requirements  of  smoothness,  then 
the  full  generality  of  shapes  cannot  be  defined  with  a  finite  number  of  parameters,  and  one  must  use  the 
concept  of  the  Ftechet  derivative  of  the  cost  with  respect  to  a  function.  Qearly,  such  a  derivative  cannot  be 
determined  directly  by  finite  differences  of  the  design  parameters  because  there  are  now  an  infinite  number  of 
these.  Using  techniques  of  control  theory,  however,  the  gradient  can  be  determined  indirectly  by  solving  an 
adjoint  equation  which  has  coefficients  defined  by  the  solution  of  the  flow  equations.  The  cost  of  solving  the 
adjoint  equation  is  comparable  to  that  of  solving  the  flow  equations.  Thus  the  gradient  can  be  determined 
with  roughly  the  computational  costs  of  two  flow  solutions,  independently  of  the  number  of  design  variables, 
which  may  be  infinite  if  the  boundary  is  regarded  as  a  free  surface.  The  underlying  concepts  are  clarified  by 
the  following  abstract  description  of  the  adjoint  method.  For  flow  about  an  airfoil  or  wing,  the  aerodynamic 


36 


Fig.  17.  Pressure  Contours  in  the  Bow  Region 


properties  which  define  the  cost  function  are  functions  of  the  flow-field  variables  (w)  and  the  physical  location 
of  the  boundary,  which  may  be  represented  by  the  function  r,  say.  Then 

I  =  I{w,r), 


and  a  change  in  T  results  in  a  change 


in  the  cost  function.  Here,  the  subscripts  7  and  JT  are  used  to  distinguish  the  contributions  due  to  the 
variation  bw  in  the  flow  solution  from  the  change  associated  directly  with  the  modification  dT  in  the  shape. 
This  notation  assists  in  grouping  the  numerous  terms  that  arise  during  the  derivation  of  the  full  Navier- 
Stokes  adjoint  operator,  outlined  in  the  next  section,  so  that  the  basic  structure  of  the  approach  as  it  is 
sketched  in  the  present  section  can  easily  be  recognized. 

Suppose  that  the  governing  equation  K  which  expresses  the  dependence  of  w  and  J-  within  the  flowfield 
domain  D  can  be  written  as 


Then  8w  is  determined  from  the  equation 


R{w,T)  =  0. 


,  \dR 
[dT 


Since  the  variation  8R  is  zero,  it  can  be  multiplied  by  a  Lagrange  Multiplier  V*  and  subtracted  from  the 
variation  81  without  changing  the  result.  Thus  equation  (43)  can  be  replaced  by 


Choosing  to  satisfy  the  adjoint  equation 
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the  first  term  is  eliminated,  and  we  find  that 


(48) 


SI  =  gsT, 


where 

The  advantage  is  that  (48)  is  independent  of  Sw,  with  the  result  that  the  gradient  of  7  with  respect  to  an 
arbitrary  number  of  design  variables  can  be  determined  withou;the  need  for  additional  flow-field  evaluations. 
In  the  case  that  (44)  is  a  partial  differential  equation,  the  adjoint  equation  (47)  is  also  a  partial  differential 
equation  and  determination  of  the  appropriate  boundary  conditions  requires  careful  mathematical  treatment. 

In  reference  [82]  the  author  derived  the  adjoint  equations  for  transonic  flows  modelled  by  both  the 
potential  flow  equation  and  the  Euler  equations.  The  theory  was  developed  in  terms  of  partial  differential 
equations,  leading  to  an  adjoint  partial  differential  equation.  In  order  to  obtain  numerical  solutions  both 
the  flow  and  the  adjoint  equations  must  be  discretized.  The  control  theory  might  be  applied  directly  to  the 
discrete  flow  equations  which  result  from  the  numerical  approximation  of  the  flow  equations  by  finite  element, 
finite  volume  or  finite  difference  procedures.  This  leads  directly  to  a  set  of  discrete  adjoint  equations  with  a 
matrix  which  is  the  transpose  of  the  Jacobian  matrix  of  the  full  set  of  discrete  nonlinear  flow  equations.  On 
a  three-dimensional  mesh  with  indices  i,j,k  the  individual  adjoint  equations  may  be  derived  by  collecting 
together  all  the  terms  multiplied  by  the  variation  of  the  discrete  flow  variable  Wjj-.fc.  The  resulting 

discrete  adjoint  equations  represent  a  possible  discretization  of  the  adjoint  partial  differential  equation.  If 
these  equations  are  solved  exactly  they  can  provide  an  exact  gradient  of  the  inexact  cost  function  which 
results  from  the  discretization  of  the  flow  equations.  The  discrete  adjoint  equations  derived  directly  from  the 
discrete  flow  equations  become  very  complicated  when  the  flow  equations  are  discretized  with  higher  order 
upwind  biased  schemes  using  flux  limiters.  On  the  other  hand  any  consistent  discretization  of  the  adjoint 
partial  differential  equation  will  yield  the  exact  gradient  in  the  limit  as  the  mesh  is  refined.  The  trade-off 
between  the  complexity  of  the  adjoint  discretization,  the  accuracy  of  the  resulting  estimate  of  the  gradient, 
and  its  impact  on  the  computational  cost  to  approach  an  optimum  solution  is  a  subject  of  ongoing  research. 

The  true  optimum  shape  belongs  to  an  infinitely  dimensional  space  of  design  parameters.  One  motivation 
for  developing  the  theory  for  the  partial  differential  equations  of  the  flow  is  to  provide  an  indication  in 
principle  of  how  such  a  solution  could  be  approached  if  sufficient  computational  resources  were  available. 
Another  motivation  is  that  it  highlights  the  possibility  of  generating  ill  posed  formulations  of  the  problem. 
For  example,  if  one  attempts  to  calculate  the  sensitivity  of  the  pressure  at  a  particular  location  to  changes 
in  the  boundary  shape,  there  is  the  possibility  that  a  shape  modification  could  cause  a  shock  wave  to  pass 
over  that  location.  Then  the  sensitivity  could  become  unbounded.  The  movement  of  the  shock,  however, 
is  continuous  as  the  shape  changes.  Therefore  a  quantity  such  as  the  drag  coefficient,  which  is  determined 
by  integrating  the  pressure  over  the  surface,  also  depends  continuously  on  the  shape.  The  adjoint  equation 
allows  the  sensitivity  of  the  drag  coefficient  to  be  determined  without  the  explicit  evaluation  of  pressure 
sensitivities  which  would  be  ill  posed. 

The  discrete  adjoint  equations,  whether  they  are  derived  directly  or  by  discretization  of  the  adjoint 
partial  differential  equation,  are  linear.  Therefore  they  could  be  solved  by  direct  numerical  inversion.  In 
three-dimensional  problems  on  a  mesh  with,  say,  n  intervals  in  each  coordinate  direction,  the  number  of 
unknowns  is  proportional  to  n®  and  the  bandwidth  to  n^.  The  complexity  of  direct  inversion  is  proportional 
to  the  number  of  unknowns  multiplied  by  the  square  of  the  bandwidth,  resulting  in  a  complexity  proportional 
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to  -nJ.  The  cost  of  direct  inversion  can  thus  become  prohibitive  as  the  mesh  is  refined,  and  it  becomes  more 
efficient  to  use  iterative  solution  methods.  Moreover,  because  of  the  similarity  of  the  adjoint  equations  to 
the  fiow  equations,  the  same  iterative  methods  which  have  been  proved  to  be  efficient  for  the  solution  of  the 
flow  equations  are  efficient  for  the  solution  of  the  adjoint  equations. 

The  control  theory  formulation  for  optimal  aerodynamic  design  has  proved  effective  in  a  variety  of 
applications  [83,  74,  159].  Pironneau  has  studied  the  use  of  control  theory  for  optimal  shape  dasign  of 
systems  governed  by  elliptic  equations  [151],  and  more  recently  the  Navier-Stokes  equations,  and  also  wave 
reflection  problems  [150].  The  adjoint  equations  have  also  been  used  by  Baysal  and  Eleshaky  [20],  and 
by  Ta’asan,  Kuruvila  and  Salas  [187],  who  have  implemented  a  one  shot  approach  in  which  the  constraint 
represented  by  the  flow  equations  is  only  required  to  be  satisfied  by  the  final  converged  solution.  In  their 
work,  computational  costs  are  also  reduced  by  applying  multigrid  techniques  to  the  geometry  modifications  as 
well  as  the  solution  of  the  flow  and  adjoint  equations.  Adjoint  methods  have  been  applied  to  incompressible 
viscous  flow  problems  by  Cabuk  and  Modi  [31,  72]  and  by  Desai  and  Ito  [49].  Recent  applications  of  adjoint 
methods  on  unstructured  meshes  include  the  work  of  Anderson  and  Venkatakrishnan  [10],  and  Elliot  and 
Peraire  [53]. 

5.3.  Three-Dimensional  Design  using  the  Compressible  Navier-Stokes  Equations.  In  order 
to  illustrate  the  application  of  control  theory  to  aerodynamic  design  problems,  this  section  treats  three- 
dimensional  wing  design  using  the  compressible  Navier  Stokes  equations  to  model  the  flow.  In  comparison 
with  incompressible  viscous  flow,  the  adjoint  equations  contain  numerous  extra  terms  which  arise  from  the 
variation  of  the  energy  equation.  In  order  to  simplify  the  calculation  of  the  effect  of  shape  changes  it  is 
convenient  to  introduce  a  body-fitted  coordinate  system,  so  that  the  flow  and  adjoint  equations  are  solved 
in  a  fixed  computational  domain.  Suppose  that  the  transformation  to  computational  coordinates  (^i,  6,  ^3) 
is  defined  by  the  metrics 

Ki^  =  {K),  = 

The  Navier-Stokes  equations  (1-3)  can  then  be  written  in  the  computational  domain  as 

(49)  ^  ^ 
where 

(50)  R(ro)  =  ^  (Fi  -  F„ J 

where  the  inviscid  and  viscous  flux  contributions  are  now  defined  with  respect  to  the  computational  cell 
faces  by  F*  =  Sijfj  and  F„i  =  and  the  quantity  Sij  =  is  used  to  represent  the  projection  of  the 

cell  face  along  the  Xj  axis.  For  convenience,  the  coordinates  describing  the  fixed  computational  domain 
are  chosen  so  that  each  boundary  conforms  to  a  constant  value  of  one  of  these  coordinates.  Variations  in 
the  shape  then  result  in  corresponding  variations  in  the  mapping  derivatives  defined  by  Kij. 

Suppose  that  the  performance  is  measured  by  a  cost  function 

1=  f  M{w,S)dBi+  f  V{w,S)dD^, 

Jb  Jv 

containing  both  boundary  and  field  contributions  where  dB^  and  dD^  are  the  surface  and  volume  elements  in 
the  computational  domain.  In  general,  MandV  will  depend  on  both  the  flow  variables  w  and  the  metrics 


dxj 
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S  defining  the  computational  space.  The  design  problem  is  now  treated  as  a  control  problem  where  the 
boundary  shape  represents  the  control  function,  which  is  chosen  to  minimize  I  subject  to  the  constraints 
defined  by  the  flow  equations  (49).  A  shape  change  produces  a  variation  in  the  flow  solution  Sw  and  the 
metrics  SS  which  in  turn  produce  a  variation  in  the  cost  function 

(51)  51  =  [  SM{w,S)dB^+  j  5V{w,S)dD^. 

JB  Jv 


Here 


5M.  =  [Adu)]/  +  5Ma, 

(52)  5V  —  \Pw\i  +  STn, 

where  we  continue  to  use  the  subscripts  I  and  II  to  distinguish  between  the  contributions  associated  with 
the  variation  of  the  flow  solution  6w  and  those  associated  with  the  metric  variations  SS.  Thus  [Mw]i  and 
[Pa,]  I  represent  ^  and  with  the  metrics  fixed,  while  6Mii  and  SVu  represent  the  contribution  of  the 
metric  variations  SS  to  SAd  and  SV. 

In  the  steady  state,  the  constraint  equation  (53)  specifies  the  variation  of  the  state  vector  Sw  by 

(53)  4-'^(Fi-F«)  =  0. 

Here  6Fi  and  5Fvi  can  also  be  split  into  contributions  associated  with  5w  and  5S  using  the  notation 


(54) 


5Fi  —  [Fi^]j  dw  +  SFiij 
SFvi  =  [Fviwli  4" 


The  inviscid  contributions  are  easily  evaluated  as 

=  5Fai„  =  5Sijfj. 

The  details  of  the  viscous  contributions  are  complicated  by  the  additional  level  of  derivatives  in  the  stress 
and  heat  flux  terms. 

Multiplying  by  a  co-state  vector  -0,  which  will  play  an  analogous  role  to  the  Lagrange  multiplier  intro¬ 
duced  in  equation  (46),  and  integrating  over  the  domain  produces 

(55)  f  V^HFi-Fai)dVi  =  0. 

Jv 

Assuming  that  is  differentiable  this  may  be  integrated  by  parts  to  give 

(56)  /  miP^S  {Fi  -  F^i)  dBi  -  [  ^5  {Fi  -  F^)  dD^  =  0. 

Jb  Jv 

This  equation  results  directly  from  taking  the  variation  of  the  weak  form  of  the  flow  equations,  where  Y’  is 
taken  to  be  an  arbitrary  differentiable  test  function.  Since  the  left  hand  expression  equals  zero,  it  may  be 
subtracted  from  the  variation  in  the  cost  function  (51)  to  give 


—  J  [^A4  -  rii'tlP' 5  {Fi  -  F«)]  dB, 


(57) 


+ 


i 


5V+^5{Fi-F^i) 


dV^. 
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Now,  since  ip  is  an  arbitrary  differentiable  function,  it  may  be  chosen  in  such  a  way  that  61  no  longer 
depends  explicitly  on  the  variation  of  the  state  vector  5w.  The  gradient  of  the  cost  function  can  then  be 
evaluated  directly  from  the  metric  variations  without  having  to  recompute  the  variation  5w  resulting  from 
the  perturbation  of  each  design  variable. 

Comparing  equations  (52)  and  (54),  the  variation  5w  may  be  eliminated  from  (57)  by  equating  all  field 
terms  with  subscript  “I”  to  produce  a  differential  adjoint  system  governing  ip 

QiIjT 

(58)  -^[Fi^-Fviwh  +  ['Pw]i  =  0  iriD. 

Oq^i 

The  corresponding  adjoint  boundary  condition  is  produced  by  equating  the  subscript  F  boundary  terms 
in  equation  (57)  to  produce 

(59)  [Fi^  -  =  [Ma]i  on  B. 


The  remaining  terms  from  equation  (57)  then  yield  a  simplified  expression  for  the  variation  of  the  cost 
function  which  defines  the  gradient 


SI  =  [  {SMn [SFi- SFvi]ij}dB^ 

Jb 

(60)  +  J  [SFi  —  SF^i]  jf  ^ 

Since  the  contributions  due  to  the  metric  variations  need  not  have  been  integrated  by  parts  this  is  equivalent 
to 

(61)  SI  =  f  SMiidB^  H-  f  {SPn  “  dD^ 

Jb  Jv 

where  SRn  is  the  variation  in  the  residual  defined  by  equation  (5.3)  due  to  variations  in  the  metrics  with 
the  fiow  solution  fixed. 

Taking  the  transpose  of  equation  (58),  the  inviscid  adjoint  equation  may  be  written  as 

(62) 

where  the  inviscid  Jacobian  matrices  in  the  transformed  space  are  given  by 


The  derivation  of  the  viscous  adjoint  terms  is 


c  -s  -^ 


simplified  by  transforming  to  the  primitive  variables 


■uF  =  {p,Ui,U2,U3,pf' , 


because  the  viscous  stresses  depend  on  the  velocity  derivatives  while  the  heat  fluxes  can  be  expressed 
as 


d  /pN 
dxi  \pj 


Here 


_  k  _  7  /X 
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where  k  is  the  conductivity,  R  is  the  gas  constant,  and  Pr  is  the  Prandtl  number.  The  relationship  between 
the  conservative  and  primitive  variations  are  defined  by  the  expressions 


Sw  =  MSui,  Sw  =  M  ^Sw, 

which  make  use  of  the  transformation  matrices  M  =  ||  and  .  The  conservative  and  primitive 

adjoint  operators  L  and  L  corresponding  to  the  variations  5w  and  5w  are  then  related  by 

f  SvFl/ipdV^^  f  SnFLipdD^, 

Jv  Jv 


L  - 

where 

'  I  Ui  U2  uz 

0/^00  pui 

M'^  =  0  0  p  0  pU2 

0  0  0  p  pm 

0  0  0  0  ^  _ 

The  detaUs  of  the  derivation  of  the  viscous  adjoint  operator  are  provided  in  [92]  with  the  simplification 
that  variations  in  the  transport  coefficients  are  ignored.  It  is  convenient  to  introduce  the  notation 

i’i+i  =(f>i,  i=  1,2,3, 'Ip5  =  0- 

in  order  to  allow  the  use  of  the  summation  convention  for  repeated  indices  over  the  range  1  to  3.  Then, 
collecting  together  the  contributions  from  the  momentum  and  energy  equations,  the  viscous  adjoint  operator 
in  primitive  variables  can  finally  be  expressed  as 


(L^)i  = 


-^—±(Srk^ 

d^i  I**"'  r  dxi. 


d  (r,  \  f  de 


for  *=1.2,3 


+  XSijUk  - 


The  conservative  viscous  adjoint  operator  may  then  be  obtained  by  the  transformation 

L  = 


The  details  of  the  formula  for  the  gradient  depend  on  the  way  in  which  the  boundary  shape  is  param¬ 
eterized  as  a  function  of  the  design  variables,  and  the  way  in  which  the  mesh  is  deformed  as  the  boundary 
is  modified.  Using  the  relationship  between  the  mesh  deformation  and  the  surface  modification,  the  field 
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integral  is  reduced  to  a  surface  integral  by  integrating  along  the  coordinate  lines  emanating  from  the  surface. 
Thus  the  expression  (26)  for  SI  is  finally  reduced  to  the  form  of  equation  (48) 


SI=  [  gSJ^dB^ 

JB 

where  T  represents  the  design  variables,  and  Q  is  the  gradient,  which  is  a  function  defined  over  the  boundary 
surface. 

The  boundary  conditions  satisfied  by  the  flow  equations  restrict  the  form  of  the  left  hand  side  of  the 
adjoint  boundary  condition  (59).  Consequently,  the  boundary  contribution  to  the  cost  fimction  M  cannot 
be  specified  arbitrarily.  Instead,  it  must  be  chosen  from  the  class  of  functions  which  allow  cancellation  of 
all  terms  containing  6w  in  the  boundary  integral  of  equation  (57).  On  the  other  hand,  there  is  no  such 
restriction  on  the  specification  of  the  field  contribution  to  the  cost  function  V,  since  these  terms  may  always 
be  absorbed  into  the  adjoint  field  equation  (58)  as  source  terms. 

The  costate  solution  -0  is  a  legitimate  test  function  for  the  weak  form  of  the  flow  equations  only  if  it 
is  differentiable.  Smoothness  should  also  be  preserved  in  the  redesigned  shape.  It  is  therefore  crucially 
important  to  introduce  appropriate  smoothing  procedures.  In  order  to  avoid  discontinuities  in  the  adjoint 
boundary  condition  which  would  be  caused  by  the  appearance  of  shock  waves,  the  cost  function  for  the 
target  pressure  pd  may  be  modified  to  the  form 


J 

XiZ 


d  ,  dz 


and  the  smooth  quantity  Z  replaces  p-pd  in  the  adjoint  boundary  condition. 

Independent  movement  of  the  boundary  mesh  points  could  produce  discontinuities  in  the  designed  shape. 
In  order  to  prevent  this  the  gradient  may  be  also  smoothed.  Both  explicit  and  implicit  smoothing  procedures 
are  useful.  Suppose  that  the  movement  of  the  surface  mesh  points  were  defined  by  local  B-splines.  In  the 
case  of  a  uniform  one-dimensional  mesh,  a  B-sphne  with  a  displacement  d  centered  at  the  mesh  point  i  would 
produce  displacements  d/4  at  *  -f  1  and  i-1  and  zero  elsewhere,  while  preserving  continuity  of  the  first  and 
second  derivatives.  Thus  we  can  suppose  that  the  discrete  surface  displacement  has  the  form 


SS  =  Bd, 


where  B  is  a  matrix  with  coefficients  defined  by  the  B-splines,  and  dk  is  the  displacement  associated  with 
the  B-spline  centered  at  i.  Then,  using  the  discrete  formulas,  to  first  order  the  change  in  the  cost  is 

SI  =  G'^SS  =  G'^Bd. 
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Thus  the  gradient  with  respect  to  the  B-spline  coefficients  is  obtained  by  multiplying  Qhy  and  a  descent 
step  is  defined  by  setting 

d  =  SS  =  Bd=  ^XBB'^g 


where  A  is  sufficiently  small  and  positive.  The  coefficients  of  B  can  be  renormalized  to  produce  unit  row 
sums.  With  a  uniform  mesh  spacing  in  the  computational  domain  this  formula  is  equivalent  to  the  use  of  a 
gradient  modified  by  two  passes  of  the  explicit  smoothing  procedure 


Gi,h  =  g0i-l,/c  + 


with  a  similar  smoothing  procedure  in  the  k  discretization. 

Implicit  smoothing  may  also  be  used.  The  smoothing  equation 


“  Gi,k)  +  ei_:L^]^{gi,k  —  Qi-l,k)  —  Gi,k 

approximates  the  differential  equation 

^  dCdi  ^ 

If  one  sets  =  —XQ,  then  to  first  order  the  change  in  the  cost  is 

gSSd^dr) 


-II 

=  -^Il(^  +  ‘  (^)  j  ■*£<'’' 


<0, 


assuring  an  improvement  if  A  is  sufficiently  small  and  positive,  unless  the  process  has  already  reached  a 
stationary  point  at  which  g  =  0. 

5.4.  Examples  of  Viscous  Design.  Due  to  the  high  computational  cost  of  viscous  design,  a  two-stage 
design  strategy  has  been  adopted  in  practical  design  applications.  In  the  first  stage,  a  design  calculation  is 
performed  with  the  Euler  equations  to  minimize  the  drag  at  a  given  lift  coefficient  by  modifying  the  wing 
sections  with  a  fixed  planform.  In  the  second  stage,  the  pressure  distribution  of  the  Euler  solution  is  used 
as  the  target  pressmre  distribution  for  inverse  design  with  the  Navier-Stokes  equations.  Comparatively  small 
modifications  are  required  in  the  second  stage,  so  that  it  can  be  accomplished  with  a  small  number  of  design 
cycles. 

Figures  25  and  26  illustrate  the  use  of  this  strategy  for  the  re-design  of  a  wing  representative  of  wide 
body  transport  aircraft.  The  design  point  was  taken  as  a  lift  coefficient  of  .55  at  a  Mach  number  of  .83. 
Figure  25  illustrates  the  Euler  redesign,  which  was  performed  on  a  mesh  with  192  x  32  x  48  cells,  displaying 
both  the  geometry  and  the  upper  surface  pressure  distribution,  with  negative  Cp  upwards.  The  initial  wing 
shows  a  moderately  strong  shock  wave  across  most  of  the  top  surface,  as  can  be  seen  in  Figure  25a.  Sixty 
dpsigu  cycles  were  needed  to  produce  the  shock  free  wing  shown  in  Figure  25b,  with  an  indicated  drag 
reduction  of  15  counts  from  .0196  to  .0181.  Figure  26  shows  the  viscous  redesign  at  a  Reynolds  number  of 
12  million.  This  was  performed  on  a  mesh  with  192  x  64  x  48  cells,  with  32  intervals  normal  to  the  wing 
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concentrated  inside  the  boundary  layer  region.  In  Figure  26a  it  can  be  seen  that  the  Enler  design  produces 
a  weak  shock  due  to  the  displacement  effects  of  the  boundary  layer.  Ten  design  cycles  were  needed  to  recover 
the  shock  free  wing  shown  in  Figure  26b.  It  is  interesting  that  the  wing  section  modifications  between  the 
initial  wing  of  Figure  25a  and  the  final  wing  of  Figure  26b  are  remarkably  small. 

These  results  were  sufficiently  promising  that  it  was  decided  by  McDonnell  Douglas  to  evaluate  the 
method  for  industrial  use,  and  it  was  used  to  support  design  studies  for  the  MDXX  project.  Reference  [88] 
provides  a  more  detailed  discussion  of  the  results  of  this  experience.  Early  in  the  project  it  became  apparent 
that  the  fuselage  effects  are  too  large  to  be  ignored  and  that  optimization  at  a  single  design  point  could  lead  to 
unsatisfactory  off-design  performance.  Therefore  we  focused  on  the  optimization  of  wing-body  combinations 
at  three  design  points.  In  viscous  design  it  was  also  found  that  there  were  discrepancies  between  the  results 
of  the  design  calculations,  which  were  initially  performed  on  a  relatively  coarse  grid  with  192  x  64  x  48  cells, 
and  the  results  of  subsequent  analysis  calculations  performed  on  finer  meshes  to  verify  the  design. 

In  order  to  allow  the  use  of  finer  meshes  with  overnight  turnaround,  the  code  was  therefore  modified  for 
parallel  computation.  Using  the  parallel  implementation,  viscous  design  calculations  have  been  performed 
on  meshes  with  1.8  million  mesh  points.  Starting  from  a  preliminary  inviscid  design,  20  design  cycles  are 
usually  sufficient  for  a  viscous  re-design  in  inverse  mode,  with  the  smoothed  inviscid  results  providing  the 
target  pressure.  Such  a  calculation  can  be  completed  in  about  l\  hours  using  48  processors  of  an  IBM  SP2. 

As  an  illustration  of  the  results  that  could  be  obtained,  Figure  27  shows  a  wing-body  design  with 
sweep  back  of  about  38  degrees  at  the  1/4  chord.  Starting  from  the  result  of  an  Euler  design,  the  viscous 
optimization  produced  an  essentially  shock  free  wing  at  a  cruise  design  point  of  Mach  .86,  with  a  lift  coefficient 
of  .6  for  the  wing  body  combination  at  a  Reynolds  number  of  101  milUon  based  on  the  root  chord.  Figure 
27  shows  the  design  point,  while  the  evolution  of  the  design  is  shown  in  Figure  28,  using  software  provided 
by  J.  Vassberg  of  Douglas  Aircraft.  In  this  case  the  pressure  contours  are  for  the  final  design.  This  wing 
is  quite  thick,  with  a  thickness  to  chord  ratio  of  more  than  14  percent  at  the  root  and  9  percent  at  the  tip. 
The  design  offers  excellent  performance  at  the  nominal  cruise  point.  Figures  29  and  30  show  the  results  of 
a  Mach  number  sweep  to  determine  the  drag  rise.  It  can  be  seen  that  a  double  shock  pattern  forms  below 
the  design  point,  while  there  is  actually  a  slight  increase  in  the  drag  coefficient  of  about  1  |  counts  at  Mach 
.85.  The  tendency  to  produce  double  shocks  below  the  design  point  is  typical  of  supercritical  wings.  This 
wing  has  a  low  drag  coefficient,  however,  over  a  wide  range  of  conditions. 

6.  Outlook  and  Conclusions.  Better  algorithms  and  better  computer  hardware  have  contributed 
about  equally  to  the  progress  of  computational  science  in  the  last  two  decades.  In  1970  the  Control  Data 
6600  represented  the  state  of  the  art  in  computer  hardware  with  a  speed  of  about  10®  operations  per  second 
(one  megaflop),  while  in  1990  the  8  processor  Cray  YMP  offered  a  performance  of  about  10®  operations  per 
second  (one  gigafiop).  Correspondingly,  steady-state  Euler  calculations  which  required  5,000-10,000  steps 
prior  to  1980  could  be  performed  in  10-50  steps  in  1990  using  multigrid  acceleration.  With  the  advent 
of  massively  parallel  computers  it  appears  that  the  progress  of  computer  hardware  may  even  accelerate. 
Teraflop  machines  offering  further  improvement  by  a  factor  of  1,000  are  likely  to  be  available  within  a  few 
years.  Parallel  architectures  will  force  a  reappraisal  of  existing  algorithms,  and  their  effective  utilization  will 
require  the  extensive  development  of  new  parallel  software. 

In  parallel  with  the  transition  to  more  sophisticated  algorithms,  the  present  challenge  is  to  extend  the 
effective  use  of  CFD  to  more  complex  applications.  A  key  problem  is  the  treatment  of  multiple  space  and 
time  scales.  These  arise  not  only  in  turbulent  flows,  but  also  in  many  other  situations  such  as  chemically 
reacting  flows,  combustion,  flame  fronts  and  plasma  dynamics.  Another  challenge  is  presented  by  problems 
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with  moving  boundaries.  Examples  include  helicopter  rotors,  and  rotor-stator  interaction  in  turbomachinery. 
It  can  be  anticipated  that  interdisciplinary  applications  in  which  CFD  is  coupled  with  the  computational 
analysis  of  other  properties  of  the  design  will  play  an  increasingly  important  role.  These  applications  may 
include  structural,  thermal  and  electromagnetic  analysis.  Aeroelastic  problems  and  integrated  control  system 
and  aerodynamic  design  are  likely  target  areas.  The  development  of  improved  algorithms  continues  to 
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Fig.  18.  Concept  for  a  numerical  wind  tunnel. 
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Fig.  19.  Advanced  numerical  wind  tunnel 

be  important  in  providing  the  basic  building  blocks  for  numerical  simulation.  In  particular,  better  error 
estimation  procedures  must  be  developed  and  incorporated  in  the  simulation  software  to  provide  error 
control.  The  basic  simulation  software  is  only  one  of  the  needed  ingredients,  however.  The  flow  solver 
must  be  embedded  in  a  user-friendly  system  for  geometry  modeling,  output  analysis,  and  data  management 
that  will  provide  a  complete  numerical  design  environment.  These  are  the  ingredients  which  are  needed  for 
the  full  realization  of  the  concept  of  a  numerical  wind  tunnel.  Figures  18  and  19  illustrate  the  way  in  which 
a  numerical  wind  tunnel  might  evolve  from  current  techniques,  which  involve  massive  data  handling  tasks, 
to  a  fully  integrated  design  environment. 

In  the  long  run,  computational  simulation  should  become  the  principal  tool  of  the  aerodynamic  design 
process  because  of  the  flexibility  it  provides  for  the  rapid  and  comparatively  inexpensive  evaluation  of 
alternative  designs,  and  because  it  can  be  integrated  in  a  numerical  design  environment  providing  for  both 
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multi-disciplinary  analysis  and  multi-disciplinary  optimization. 


Author’s  Note.  This  paper  was  presented  as  the  second  Theodorsen  Lecture  at  ICASE.  I  warmly 
appreciate  the  invitation  to  give  this  lecture.  Theodorsen’s  solution  of  the  potential  flow  past  arbitrary 
profiles  via  conformal  mapping  was  probably  the  flrst  successful  application  of  a  numerical  method  to  applied 
aerodynamics,  and  it  provided  the  basis  for  the  development  of  the  NACA  five  digit  series  of  airfoils.  The 
second  half  of  the  lecture  addresses  the  problem  of  aerodynamic  shape  design.  One  of  the  many  contributions 
of  Lighthill,  who  gave  the  flrst  Theodorsen  lecture,  was  the  solution  of  the  inverse  problem  of  airfoil  design 
for  a  specified  pressure  distribution  via  conformal  mapping.  The  author’s  first  application  of  design  using 
control  theory  [82,  83]  extends  Lighthill’s  method  to  design  for  transonic  potential  flow. 

The  research  described  in  this  paper  has  benefited  greatly  from  the  generous  support  of  the  AFOSR 
imder  grant  number  AFOSR  91-0391,  and  ARPA  under  grant  number  N00014-92-J-1976.  I  am  grateful  to 
Juan  Alonso,  Luigi  MartenelU,  JoAnn  Love,  and  Stephen  Guattery  for  their  assistance  in  assembling  the 
text  with 
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Fig.  20.  Navier-Stokes  Predictions  for  the  F-18  Wing-Fuselage  at  Large  Incidence 


22a:  Cp  after  35  Cycles. 
Cl  =  0.3654,  Cd  =  0.0232. 


Figure  22b: 


Convergence. 


Fig.  22.  NACA-0012  Airfoil  at  Mach  0.800  andot  =  1.25° H-CUSP  Scheme. 


d 


23a:  Cp  after  35  Cycles. 
Cl  =  0.3861,  Cd  -  0.0582. 


Figure  23b:  Convergence. 


Fig.  23.  NACA-0012  Airfoil  at  Mach  0.850  anda  =  1.0° H-CUSP  Scheme. 


24a:  12.50%  Span. 

Cl  =  0.2933,  Cd  =  0.0274. 


24c:  50.00%  Span. 

Cl  =  0.3262,  Cd  =  0.0089. 


24b:  31.25%  Span. 

Cl  =  0.3139,  Cd  =  0.0159. 


24d:  68.75%  Span. 

Cl  =  0.3195,  Cd  =  0.0026. 


Fig.  24.  Oneva,  M6  W^ing,  Mach  0.840,  Angle  of  Atta,ckS.0G>^ ,  192x32x48  Mesh.  Cl  —  0.3041,  Co  —  0.0131.  H-CUSP 
scheme. 


62 


Initial  Wing.  Cp  on  Upper  Surface. 

Figure  25a:  M  =  .83,  Cj  =  .5498,  Cd  =  .0196,0:  =  2.410°. 


Redesigned  wing.  Cp  on  Upper  Surface. 

Figure  25b:  M  =  .83,  Ci  =  .5500,  Cd  =  .0181,  o  =  1,959°. 

Fig.  25.  Redesign  of  the  wing  of  a  wide  transport  aircraft.  Stage  1  Inviscid  design  :  60  design  cycles  in  drag  reduction 
mode  with  forced  lift. 


Initial  Wing.  Cp  on  Upper  Surface. 

Figure  26a:  M  =  0.83, Q  =  .5506,  Cd  =  .0199,0;  =  2.317° 


Redesigned  wing.  Cp  on  Upper  Surface. 

Figure  26b:  M  =  0.83,  Ci  =  .5508,  Ca  =  .0194,  o  =  2.355° 

Fig.  26.  Redesign  of  the  wing  of  a  wide  transport  aircraft.  Stage  2:  Viscous  re-design.  10  design  cycles  in  inverse  mode 
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COMPARISON  OF  CHORDWISE  PRESSURE  DISTRIBUTIONS 
MPX5X  WING-BODY 


Fig.  27.  Pressure  distribution  of  the  MPX5X  at  its  design  point. 


COMPARISON  OF  CHORDWISE  PRESSURE  DISTRIBUTIONS 
MPX5X  WING-BODY 


Fig.  28.  Optimization  Sequence  in  the  design  of  the  MPX5X. 
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COMPARISON  OF  CHORDWISE  PRESSURE  DISTRIBUTIONS 
MPX5X  WING-BODY 


Fig.  29.  Off  design  performance  of  the  MPX5X  below  the  design  point. 


COMPARISON  OF  CHORDWISE  PRESSURE  DISTRIBUTIONS 
MPX5X  WING-BODY 


Fig.  30.  Off  design  performance  of  the  MPX5X  above  the  design  point. 
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